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I.  INTRODUCTION 


A  finite  element  model,  (FEM)  is  a  computational  representation  of  a  structure  or 
object  that  is  generally  used  to  solve  a  complex  problem  for  which  a  closed  fonn  solution 
is  unavailable.  The  nature  of  the  problem  frequently  means  that  these  structures  or 
objects  are  very  complex  in  design.  For  many  systems,  the  FEM  is  the  only  means  by 
which  the  response  of  the  system  can  be  calculated. 

The  FEM  requires  a  number  of  parameters,  including,  but  not  limited  to,  physical 
geometry,  density  and  modulus  of  elasticity.  Although  a  structure  is  defined  by  a  large 
number  of  parameters  only  a  small  number  of  parameters  can  be  measured  in  a  modal 
test.  It  is  these  measured  parameters  from  a  modal  test  are  the  modal  parameters  of  a 
structure,  i.e.,  natural  frequency  and  mode  shapes.  It  is  these  parameters  that  define  a 
structure  and  will  be  the  focus  of  this  paper.  The  ability  of  a  model  to  accurately  predict 
the  modal  parameters  is  of  vital  importance.  Often  the  model  must  be  corrected  or 
updated  in  order  to  accurately  reflect  the  correct  parameters.  To  update  the  modal  real 
world  data  must  be  applied  to  and  correct  the  FEM. 

In  most  real  world  problems  a  limited  number  of  measured  modal  parameters 
from  the  structure  creates  an  underdetermined  problem.  This  underdetennined  problem 
is  a  product  of  the  large  number  of  adjustable  parameters  compared  to  the  relatively  small 
number  of  measured  parameters.  Often  the  number  of  measured  parameters  are  limited 
due  to  equipment  limitations  (i.e.,  inability  to  measure  modes  at  high  frequencies). 
However,  when  measuring  and  obtaining  data  from  a  structure,  the  number  of 
measurements  can  also  be  limited  by  the  number  of  measuring  devices  available  (i.e., 
accelerometers).  Thus,  it  is  necessary  to  place  the  devices  in  either  locations  of  interest 
or  locations  of  ease.  These  measured  locations  result  in  the  “analysis  set”  of  coordinates 
or  “ASET”.  The  remaining  set  is  referred  to  as  the  “omitted  set”  or  “OSET”.  Due  to 
most  finite  element  models  having  a  large  number  of  coordinates,  or  degrees  of  freedom 
(DOF),  the  number  of  OSET  coordinates  is  much  greater  than  the  number  of  ASET 
coordinates. 


1 


One  method  of  solving  the  above  problem  is  to  obtain  more  measured  parameters 
from  the  structure.  One  approach  to  this  is  to  obtain  modal  parameters  for  a  larger 
number  of  modes.  This  is  often  prohibitive  due  to  the  difficulties  associated  with 
measuring  higher  frequency  modes  for  a  structure.  Another  approach  involves  the 
application  of  additional/new  boundary  conditions  to  the  structure  and  obtaining  more 
data.  However,  this  is  often  costly  or  simply  impossible.  Thus,  a  more  feasible  solution 
method  is  desired.  One  option  is  the  method  of  Artificial  Boundary  Conditions. 

The  method  of  Artificial  Boundary  Conditions  applies  a  pseudo  boundary 
condition  to  the  structure  and  extracts  modal  parameters  from  the  pre-existing  data.  The 
term  pseudo  is  used  because  an  actual  physical  boundary  condition  is  not  applied,  the 
model  merely  reacts  as  if  additional  degrees-of-freedom  (DOF)  are  restrained.  This 
results  in  the  ability  to  extract  additional  modal  parameters  without  having  to  conduct 
additional  testing  on  the  structure.  When  applying  this  method  to  model  updating  it  can 
often  help  to  resolve  an  underdetermined  problem. 

Sensitivity  based  updating  makes  use  of  the  modal  parameters  to  help  locate  a 
difference  between  a  FEM  and  an  existing  structure.  This  difference  could  exist  for  a 
variety  of  reasons.  The  difference  could  be  due  to  damage  or  an  unknown  alteration  of 
the  structural  system  that  is  desired  to  be  located.  The  difference  could  also  be  the  result 
of  a  poorly  representative  model,  or  even  the  result  of  a  model  that  is  not  capable 
replicating  the  physics  of  the  structure.  Trying  to  locate  this  difference  is  a  problem  that 
is  generally  underdetermined  due  to  the  number  of  parameters  that  can  be  altered  and  the 
limited  number  of  measured  parameters.  With  the  application  of  ABC  a  better  solution  to 
this  problem  can  often  be  formulated. 
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II.  THEORY 


The  general  equation  of  motion  will  be  the  fundamental  starting  point  for  the 
problem. 

[M]Jij  +  [C]hJ+[X]{i}  =  {F(0}  (1) 

This  problem  consists  of  mass,  damping  and  stiffness  parameters,  where  [if], 
[M]  and  [C]  are  the  symmetric  stiffness,  mass,  and  damping  matrices,  of  size  n  x  n,  and 
|F(Y)}  is  the  excitation  vector  of  size  nxl.  They  can  be  more  explicitly  stated  in  tenns  of 
their  components  within  each  matrix,  and  are  so  stated  below. 


The  above  equation  is  for  n  number  of  DOF.  As  mention  previously,  in  real 
world  analysis,  it  is  improbable  that  all  DOF  will  be  instrumented  and  recorded.  In  this 
testing  the  group  of  DOF  that  are  instrumented  and  recorded  are  considered  the  analytical 
set  or  “ASET”  while  the  remaining  DOF  are  considered  the  omitted  or  “OSET”. 

A.  ANALYSIS  AND  OMITTED  COORDINATE  SETS 

When  conducting  modal  testing  the  system  being  analyzed  can  be  broken  up  into 
two  types  of  coordinate  sets.  Starting  with  the  EOM  listed  above  and  applying  a  simple 
hannonic  excitation 

{F(0}  =  |f|^  (3) 

Giving  a  steady  state  response  of 

{x(0}  =  jxje^  (4) 
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Now  substituting  this  into  the  EOM  for  the  n  DOF  system. 

"IX  -  r mu  -  Mj  i  rcn  -  cji  r^(o] 

:  \  :  -Q2  :  •.  :  +  y'Q  :  •.  :  +<  >  =  <  :  >  (5) 

U„i  -  [Mnl  •••  M,  J  [cHl  -  C,Jj  [xH(t)\  U, 

Taking  this  example  and  now  thinking  if  it  in  tenns  of  ASET  and  OSET  coordinate 

systems  we  can  rewrite  the  [/< ] ,  [M] ,  and  [C]  in  those  terms,  leading  to 


Looking  at  the  above,  there  obviously  exists  a  relationship  between  the  ASET  and 
OSET  coordinates  in  which  the  ASET  coordinates  can  be  used  to  develop  the  solution  for 
the  OSET  coordinates  (Gordis,  1993).  Ignoring  damping  and  assuming  that  there  is  no 
excitation  acting  on  the  OSET  coordinate  system  then  equation  can  be  written  as 


The  matrices  can  easily  be  broken  into  the  two  separate  equations.  Looking  at  the 
two  equations,  the  solution  can  be  determined  for  the  OSET  coordinates  in  terms  of 
ASET  coordinates  using  just  one  of  the  two. 

K.  Ik  ( + Jk } -  n2  lMm  ]{*, } + [m„  ]k }] = o  (8) 

Then  grouping  the  similar  terms  and  solving  in  terms  of  {xG},  Eqn.  (8)  becomes 

k  ( =  [i -n2K;'oMoa )'  [-  k£k„  +  n2K  \x„ }  (9) 

In  Eqn.  (9)  the  [i  -Q2K^Moo  j  term  can  be  written 
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where  Det[...]  indicates  the  determinate  and  Adj[...]  indicates  the  adjoint  matrix.  From 
Eqn.  (10)  it  is  evident  that  the  solution  for  {x0}  in  Eqn.  (9)  does  not  exist  at  those 
frequencies,  which  satisfy 

Det[l-n-K-J:Mj=0  (11) 

The  frequencies  which  satisfy  Eqn.  (11)  are  the  eigenvalues  of  the  system  defined 
by  the  matrices  [ K00\  and  [ M00\ .  Since  both  of  these  matrices  are  composed  solely  of 
OSET  coordinates,  the  resulting  system  is  considered  the  OSET  system  and  the  resulting 
frequencies  are  the  OSET  frequencies.  Knowing  that  the  eigenvector  solutions  to  a 
system  are  defined  by  the  unrestrained  DOF  it  should  be  noticed  that  the  OSET  system  is 
considered  the  unrestrained  DOF  and  the  ASET  the  restrained  or  grounded  set. 

The  use  of  a  limited  number  of  analyzed  DOF  to  represent  a  complete  structure  is 
imperative  in  structural  dynamics.  Even  if  the  number  of  ASET  coordinates  is  relatively 
large  it  would  still  be  far  less  than  the  literally  infinite  OSET  coordinates  (Gordis,  1999). 
Thus  even  though  a  modal  test  is  limited  in  the  number  of  ASET  coordinates,  it  is 
imperative  that  the  use  of  this  incomplete  coordinate/data  set  be  used  to  represent  a  full 
system. 

B.  SPATIALLY  INCOMPLETE  DATA 

As  mentioned  previously  a  finite  element  model  can  represent  a  real  world  system 
with  an  infinite  number  of  DOF.  However,  real  world  limitations  allow  only  a  limited 
number  of  DOF  to  be  measured  in  a  modal  test  of  a  system.  This  lack  of  a  complete  data 
for  all  DOF  creates  a  spatially  incomplete  problem.  This  leads  to  the  necessity  of  finding 
a  solution  using  the  spatially  incomplete  data  set.  One  method  to  do  this  is  using  the 
spatially  incomplete  data  to  solve  for  the  OSET  frequencies  which  will  be  derived  in  the 
next  section.  An  understanding  of  how  this  spatially  incomplete  data  compares  to 
spatially  complete  data  follows. 

Consider  first  the  complete  FRF  matrix,  which  is  n  x  n  and  then  suppose  that  the 
description  of  the  system  is  limited  to  only  certain  measured  coordinates,  thus  ignoring 
what  happens  at  the  other  coordinates.  (Note  that  ignoring  coordinates  is  not  the  same  as 
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supposing  that  the  other  coordinates  do  not  exist.)  The  resulting  model  is  now  of  the 
order  N  x  N.  It  is  clear  that  because  the  basic  system  has  not  been  altered  and  that  it  still 
has  the  same  number  of  DOF,  even  though  it  was  decided  not  to  describe  all  of  the  DOF, 
the  elements  which  remain  in  the  reduced  FRF  matrix  are  identical  to  the  corresponding 
elements  in  the  full  n  x  n  matrix.  (Ewins,  1982)  Essentially  this  is  saying  that  the 
reduced  matrix  is  composed  of  elements  of  interest  from  the  original  matrix  and  leaving 
behind  those  that  are  ignored.  This  matrix  however  still  contains  all  the  modal 
information  of  a  fully  described  matrix. 

Specifically  looking  at  the  before  mentioned  description  of  a  reduced  matrix  and 
that  the  data  contained  in  a  measured  FRF  matrix,  which  is  composed  of  response 
information  from  a  limited  number  of  measurements,  it  can  be  seen  that  it  contains  all  of 
the  infonnation  of  the  entire  system.  Also  given  that  the  measured  FRF  matrix  implicitly 
defines  a  dynamically  reduced  impedance  model  the  reduction  of  the  FE  model  is 
pursued  is  the  same  manner  (Gordis,  1996).  Therefore  the  use  of  a  spatially  incomplete 
data  set  to  solve  for  the  entire  system  is  justified.  This  is  of  particular  significance  in  the 
application  of  FRFs. 
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III.  FREQUENCY  RESPONSE  FUNCTION 


A  Frequency  Response  Function,  or  FRF,  is  the  ratio  of  the  output  response  of  a 
structure  due  to  an  applied  force.  The  FRF  contains  the  amplitude  and  phase  for  the 
response  for  a  specific  DOF  due  to  an  input  force  of  unit  amplitude  at  a  specific  DOF. 
This  amplitude  and  phase  describe  where  the  natural  frequencies  of  a  system  lie  and  are 
of  great  use  in  model  updating.  Additionally  FRF  data  can  be  estimated  from  output  only 
measurements  on  operating  machinery  (Hanson,  2005).  Due  to  these  attributes  the  use 
of  FRF  in  structural  dynamic  modal  testing  is  widespread. 

An  FRF  is  constructed  by  measuring  both  the  applied  force  and  the  response  of 
the  structure  to  the  applied  force  simultaneously.  The  measured  data  are  then 
transfonned  from  the  time  domain  to  the  frequency  domain  using  Fast  Fourier  Transfonn 
(FFT)  algorithms.  Due  to  this  transfonnation  the  function  ends  up  being  constructed  in 
tenns  of  real  and  imaginary  components  or  magnitude  and  phase  components. 

The  FRF  matrix,  [//(Q)] ,  relates  the  amplitude  of  stead  state  harmonic  forcing  at 
DOF  “/’  to  the  amplitude  of  the  steady  state  hannonic  response  at  DOF  “/”.  The  inverse 
of  the  FRF  matrix  is  known  as  the  impedance  matrix  [Z(Q)] . 

[Z(Q)]  =  [fT(Q)]  '  (12) 

where 

[Z(Q)]  =  [if-Q2M  +  yQC]  (13) 

A.  REDUCED  ORDER  MODEL 

The  limited  number  of  instruments  recording  data  in  a  model  describes  a  reduced 
order  model.  This  reduced  order  model’s  impedance  is  non-linearly  dependent  on  the 
impedance  of  the  full  order  model  (Gordis,  1999).  A  full  impedance  matrix  of  infinite 
DOF  is  described  below 


7 


(14) 


H 

H 


aa 

oa 


H 

H 


ao 

oo 


This  is  the  combination  of  measured  coordinates,  also  considered  experimental 
set,  a  and  omitted  coordinate  set,  o  .  The  experimental  set  can  also  be  thought  of  as  the 
reduced  model  indicated  by  the  overbar  as  shown  below. 


[wj 


(15) 


From  the  previous  discussion  of  the  impedance  and  FRF  matrix  it  is  known  that 


Zaa 

^oa 

~Haa 

HoaA 

'1 

0" 

ao 

Zoo_ 

Hao 

Hoa  J 

0 

1 

(16) 


Writing  the  matrix  in  equation  form 


ZaaHaa  +  ZaoHoa  -  1 

ZH„  +  ZH^  =  0 


(17) 


Now  solving  for  \Haa  ]  gives 

(18) 

Knowing  that  [Z]  =  [ ~K  -  D.2  M  -  jD.C~^  or  [k  -Q2m]  if  damping  is  zero  and  that 
the  eigenvalue  solution  is  defined  by  the  solution  to  [k  -Q2m]=0  it  can  be  seen  that  the 
solution  to  [H aa  ]  =  \z,aa  -  ZaoZ~^Zou  ]  *  is  defined  by  the  Zj  term.  This  shows  that 
elements  of H~la  will  be  singular  at  the  OSET  natural  frequencies  (Gordis,  1999).  To 
show  the  importance  of  this  singularity  a  brief  derivation  of  the  FRF  follows. 


B.  DERIVATION  OF  FREQUENCY  RESPONSE  FUNCTION 

Looking  at  a  impedance  matrix  in  terms  of  an  n  DOF  system 


Zl\  • 

••  V 

*i(0 

: 

>  =  < 

A 

Z«1  • 

Znn  _ 

xn(t) 

fn , 

(19) 
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Then  in  terms  of  the  FRF  matrix  for  the  same  n  DOF  system 


x,  (t) 

~Hn  • 

1 

R 

af 

7i' 

x,\t) 

*  — 

1 

^  ... 

a 

-  Hnn_ 

: 

(20) 


Knowing  that  {x}  =  ,  where  ®  is  the  mass  nonnalized  mode  shapes  of  the 

system  and  q  gives  the  generalized  coordinate  set 

[zM «}={/}  (21) 

Then  pre-multiplying  both  sides  by  [®]r 

[0]r[Z][®]{?}=[0]r{F}  (22) 


Then  expanding  Z ]  from  Eqn  (13). 


[®f  [K][0]  -Q2  [®f  [M][®]  +  7Q[®f  [C][0]  { q }  =  [®f  {F}  (23) 


Invoking  that  Eqn.  (23)  is  mass  nonnalized  [®]r[M][®]  =  1  and  assuming 


proportional  damping. 


C  0 
0  C 


Eqn  (23)  reduces  to 


\a>?  -Q2  =  [®f  {F}  (24) 


Where  a>i  coi  is  the  zth  natural  frequency  of  the  system,  C,  is  the  systems  damping 
ratio,  and  Q  is  the  forcing  frequency.  At  this  point  -Q2  +2  j^D.co.j  is  known  as  the 
modal  impedance  matrix  and  is  diagonal. 

Inverting  the  modal  impedance  matrix  to  find  the  modal  frequency  response 
matrix  gives 


W=M 


—  d.  +2j  CFl  ojj 


[®f{F} 


(25) 
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Transforming  Eqn.  (25)  back  into  physical  coordinates,  and  pre-multiplying  by 
[®]  and  making  use  of  {3}  =  [®]r  {F} 


=[*] 


ccr  -  Q2  +  2  j<^O.coi 


KM =[®M 


(26) 


From  Eqn.  (26)  [//(Q)]  can  be  extracted  and  defined  as 


[ff(n)H®], 


,[*]' 


\ar  -Q2  +2jCD-0)][ 

The  FRF  Matrix  can  also  be  written  in  summation  form  for  each  element. 


(27) 


mode 


®A  m 


[//(«)]=  2-. 

T=i  cok  -Q  +2 jgl(ok 

This  summation  form  can  be  descriptive  in  terms  of  ^//y(Q)J  . 


(28) 


mod  es 


KH=  z 


cok  -O  +  2jCD.cok 


(29) 


From  Eqn.  (29)  it  can  be  seen  that  the  FRF  contains  the  mode  shape  data  for  the 
system.  In  fact  it  should  be  noted  that  any  row  of  the  FRF  contains  all  the  mode  shape 
data  for  the  entire  system. 

The  FRF  output  display  in  the  frequency  domain  also  contains  and  displays  the 
resonance  and  anti-resonance  frequencies  of  the  system.  It  is  this  resonance  and  anti¬ 
resonance  frequency  information  that  will  be  of  interest  in  later  sections.  It  should  be 
noted  that  the  subscripts  for  the  mode  shapes  indicate  the  DOF  that  the  force  is  being 
applied  (j),  and  response  at  DOF  measured  (i).  When  these  two  are  the  same  DOF  is 
tenned  the  driving  point,  or  “point  mobility”,  function  and  has  some  unique 
characteristics.  When  ( / )  and  (/)  are  different  DOFs  it  is  termed  “transfer  mobility” 
function  (Ewins,  1982). 
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The  following  example  is  a  two  DOF  system  with  no  damping.  This  example 
comes  from  Ewins,  1982.  and  displays  the  characteristics  of  both  a  driving  point  and 
transfer  function  for  the  complete  FRF. 


Figure  1.  2  DOF  Example 


With  the  stiffness  matrix  defines  as  [if] 


Kl  +  K2  -  K2 

-k2  k2+k2 


and  the  mass  matrix  [M] 


Ml  0 

.  0  m2 


For  this  example  Ml  =  M2=  1.0  and  Kx  =  K3  =  0.4  and  K2=. 8.  This  yields  the 
following: 


Mode  shapes:  {o 1 }  = 


{-°'7071l{®n=F0'7071} 

[-0.7071J  1  J  [0.7071  J 


Natural  frequencies:  {co(rad  /  sec)} 


JO. 6325) 
Jl  .0954 J 


Applying  this  with  Eqn.  (29)  \h&  (Q)]  = 


h  ork  -Q2  +  2jCD.0)k 


For  the  driving  point  function 

\H  m)l=fefeT  ,  fefef  {-0J07j{-0J07f  |  {- 0.707 jj- 0.707 f 
11  cof-Q2  a>2  -Dr  0.4 -Q2  1.2 -Q2 

and  for  the  transfer  function. 

\„  (P-)l  bfer  ,  hf|jl’:}r  {- 0. 707 1}{- 0.707 l}r  |  i-0.707l}{0.707l}r 
L  121  O)2  -  Q2  co2-n2  0.4  -  Q2  1.2 -Q2 
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For  the  driving  point  function  it  can  be  seen  that  the  two  terms  will  have  the  same  sign 
for  forcing  frequencies  above  and  below  the  two  natural  frequencies  (resonance  peaks). 
This  creates  an  additive  effect  where  the  two  terms  are  summed  to  create  a  larger  third. 
This  is  displayed  in  the  top  portion  of  Figure  (2)  which  is  the  Driving  Point  FRF  [i/n] . 

However,  in  between  the  two  natural  frequencies  the  two  terms  are  of  opposite  sign  and 
thus  are  subtractive  in  nature.  Due  to  this  there  exists  a  forcing  frequency  that  will  cause 
the  terms  to  equal  each  other  and  thus  create  an  antiresonance.  It  should  be  noted  that  the 
“y”  axis  is  log  based  thus  the  addition  and  subtraction  of  the  two  terms  does  not  fall 
equally  between  them  on  the  plot 

Looking  at  the  transfer  function  it  is  evident  that  the  two  terms  will  have  opposite 
signs  above  and  below  the  two  natural  frequencies.  This  creates  a  subtractive  effect  in 
the  two  regions.  This  is  shown  in  the  lower  plot  of  Figure  (2)  which  is  the  H12. 
Conversely  in  the  region  between  the  two  natural  frequencies  the  two  terms  are  additive 
and  thus  do  not  create  an  antiresonance. 


Driving  Point  FRF,  H11 


Transfer  FRF,  HI 2 
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Figure  2.  2  DOF  Frequency  Response  Function 
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This  principles  demonstrated  here  can  be  extended  to  any  number  DOF  and  thus 
create  a  fundamental  rule  that  has  great  value.  That  is  that  if  two  consecutive  modes  have 
the  same  sign  for  the  modal  constants,  then  there  will  be  an  anti-resonance  frequency 
between  the  two  natural  frequencies  of  those  two  modes  (Ewins,  1982).  This  concept  of 
the  existence  of  an  antiresonance  between  any  two  modes  of  a  driving  point  FRF  is  of 
great  significance  in  the  use  of  Artificial  Boundary  Conditions. 


13 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


14 


IV.  ARTIFICIAL  BOUNDARY  CONDITIONS 


As  discussed  previously  the  ability  of  a  finite  element  model  to  accurately 
represent  the  dynamic  behavior  of  a  structure  is  the  desired  end  result.  Of  particular 
importance  are  the  similarities  between  the  natural  frequencies  of  the  model  and  the 
actual  structure.  Often  there  exists  a  disparity  between  the  two  and  it  is  desired  to  correct 
or  “update”  the  model  to  better  represent  the  structure.  In  order  to  do  this  it  is  necessary 
to  obtain  as  much  data  from  modal  testing  as  possible.  This  could  be  resolved  by 
conducting  multiple  tests  of  the  system  under  different  boundary  conditions.  This  would 
obviously  be  a  very  expensive  and  time  consuming  procedure.  The  use  of  ABC  can 
alleviate  the  need  for  additional  testing  by  mathematically  manipulating  the  existing  data. 

A.  ARTIFICIAL  BOUNDARY  CONDITIONS  DEFINED 

The  term  artificial  boundary  condition  describes  the  constraining  or  “pinning”  of 
a  specified  DOF  on  a  finite  element  model.  The  tenn  artificial  portrays  that  there  is  no 
actual  physical  boundary  conditions  applied  to  the  structure  under  test.  However,  the  test 
data  from  the  structure  is  imposed  with  computational  constraints.  The  existing  model  is 
merely  updated  with  new  boundary  conditions  that  are  easily  imposed  on  a  finite  element 
model.  This  then  generates  additional  data  that  can  be  used  in  the  comparison  of  the  two. 
Numerous  ABC  can  be  applied  to  a  single  FE  model  thus  generating  several  models 
varying  only  in  the  boundary  conditions.  The  new  boundary  condition  provides  a  new 
configuration  to  the  model,  yet  only  one  set  of  measured  test  data  is  required.  These  ABC 
are  the  boundary  conditions  that  define  the  OSET.  As  previously  shown  OSET  provide 
additional  frequency  information  about  the  model  and  this  additional  information  can 
help  to  eliminate  ill-conditioning  in  a  solution  (Gordis,  1999). 

A  simple  example  is  the  spectrum  of  antiresonance  for  any  driving-point 
frequency  response  function.  In  this  spectrum  the  antiresonance  frequencies  correspond 
to  the  mode  frequencies  of  the  structure  with  the  driving-point  DOF  constrained. 
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1.  Simple  Two  Degree  of  Freedom  Example 


A  simple  two  DOF  example  from  (Gordis,  1999)  will  demonstrate  this.  The 
system  is  composed  of  masses,  Ml  and  M2,  and  springs  K1  and  K2  and  is  undamped. 
This  system  is  shown  in  Figure  (3). 


Mi 

/\/ 

m2 

Figure  3.  Two  DOF  system 


From  Eqn.  (29)  the  driving  point  FRF  is 


modes 


k= 1 


($f)2 

o)2k -Cl2 


(32) 


Where  </>■  is  the  mass  nonnalized  mode  shape  element,  wr ,  is  the  rth  natural 
frequency,  and  Q  is  the  forcing  frequency. 

Solving  for  the  antiresonance  forcing  frequency  of  Hu(Cl)  leads  to 


Q 


2 

anti-res 


R'uco;  +  Rffi)2 

*n+*n 


(33) 


Where  the  modal  residual  is  given  by  RXJ  =  (fitff.  .  Using  values  of  Mi  =M2=1.0, 
and  Ki  =  K2=1.0.  This  generates  mass  and  stiffness  matrices 


-Kx 

“1 

-f 

-K, 

Kx+K2 

0 

2 

M 


Mi 

0  1 

“1 

0" 

0 

M2'_ 

0 

1 

yielding  natural  frequencies:  [co(rad  /  sec)} 


JO. 618) 
[1.618  J 


and  a  single  anti-resonance 


located  at  the  frequency  Qanti-res  =  V 2  rad/sec.  The  driving  point  FRF,  [HI  1]  is  shown  in 
Figure  4.  The  single  anti-resonance  is  noticeable  at  V' 2  rad/sec. 
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Driving  Point  FRF,  H11 
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Figure  4.  Hn  Driving  point  FRF. 

The  same  system  is  analyzed,  however  with  DOF  1  constrained  or  pinned  as 
shown  in  Figure  (5). 


Figure  5.  Two  DOF  system  with  ABC  employed. 

The  resulting  systems  natural  frequency  lies  at  V2  rad/sec,  which  is  identical  to 
the  original  systems  antiresonance  frequency.  This  correlation  is  shown  in  Figure  6 
below. 
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Plot  A.  Driving  Point  FRF,  H11 


40 


_40  I _ i _ i _ i _ i _ i _ 1 _ I _ i _ i _ 

0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 

rad/ sec 


t5 


o 

$ 


Plot  B.  Z  =inv(H11)  of  the  ABC  System 
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Figure  6.  Plot  A,  Driving  Point  FRF  of  system  1,  Plot  B,  Driving  Point  FRF  of 

system  2. 


2.  ABC  Frequencies  for  a  Simply  Supported  Beam 

The  following  example  from  (Fernandez,  2004)  illustrates  the  use  of  ABC  on  a 
cantilever  beam  model.  The  model  consists  of  ten  elements,  with  two  DOF  per  element 
yielding  20  total  DOF. 


Figure  7.  10  element  cantilever  beam 

Using  Eqn.  (29)  the  driving  point  FRF  was  calculated.  The  first  eight  natural 
frequencies  of  the  system  were  calculated  to  be  4.9186,  30.826,  86.332,  169.29,  280.29, 
419.91,  589.15,  789.30,  all  in  FIZ.  The  respective  anti-resonances  calculated  are  6.8697, 
43.856,  124.28,  245.53,  406.42,  586.39,  704.69,  in  Hz.  The  driving  point  FRF  of  the 
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[H33]  was  calculated  using  Eqn.  (29)  and  plotted  versus  frequency.  Next  an  ABC  was  the 
applied  to  the  model  at  DOF  3  as  shown  in  Figure  8. 


Figure  8.  10  element  cantilever  beam,  ABC  applied  at  DOF  3. 


Again  the  natural  frequencies  were  calculated,  this  time  using  the  reduced  order  [K],  [M] 
from  the  ABC  system.  The  natural  frequencies  of  the  ABC  system,  DOF  3  pinned,  under 
800  Hz  are  6.8697,  43.856,  124.28,  245.53,  406.42,  586.39,  704.69.  Comparing  these  to 
the  peaks  of  [H33]  _1  a  strong  relationship  can  be  seen.  The  same  as  example  the  plots  of 
[H33]  and  [H33]  _1  are  combined  on  Figure  (9)  for  easier  comparison  of  the  location  of 
peaks  and  anti-resonances. 


Plot  A.  Driving  point  FRF,  H33  for  Cantilever  beam 


Plot  B.  OSET  ffeq  for  ABC  system,  ASET=  [3]  pinned 


Figure  9.  Plot  A.  Driving  Point  H33  (Q)  of  free-free  beam 
Plot  B.  [H33  (Q)]'1  with  ABC  applied  at  DOF  3 
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This  demonstrates  that  the  use  of  OSET  natural  frequencies  can  be  used  in 
improving  and  underdetermined  problem  solution.  By  making  use  of  a  just  a  single  set  of 
data  from  a  modal  parameter  test  of  a  structure,  additional  data  can  be  extrapolated.  This 
additional  data  is  withdrawn  from  the  OSET  natural  frequencies  and  easily  applied  to  the 
model  updating  problem. 
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V.  SENSITIVITY  BASED  UPDATING 


As  stated  previously  the  matching  of  modal  parameters  between  a  FEM  and  the 
structure  it  represents  is  of  great  significance  in  structural  dynamics.  In  order  to  do  this 
the  model  often  must  be  updated  to  obtain  the  actual  modal  parameters  of  the  structure. 
When  updating  a  finite  element  model  to  more  closely  match  the  actual  structure  the 
method  of  sensitivity  based  updating  is  often  employed.  The  governing  equation  for 
sensitivity  based  updating  is 

|  Aw2}  =[r]{ADv}  (34) 

where  {Aw}  is  the  change  in  natural  frequency,  {ADv}  is  the  change  in  design  variable  to 
be  solved  for,  and  \T ]  is  the  sensitivity  matrix,  composed  of  first-order  sensitivities. 


A.  SENSITIVITY  MATRIX 


The  sensitivity  matrix  can  be  defined  as  how  a  change  in  a  model  parameter  for 
one  element  affects  the  other  elements  of  the  model.  This  change  could  be  a  change  of 
element  mass  or  stiffness  and  its  respective  change  on  the  systems  change  in  natural 
frequency.  In  matrix  fonn: 


bl 


8a>; 

dDv 


(35) 


where  each  column  represents  an  element  of  the  model  and  each  row  represents  a  mode 
of  the  model. 

The  programs  used  to  compose  the  sensitivity  matrix  for  this  study  did  so  doing 
the  following  steps: 

1 .  A  small  perturbation  (1%)  of  mass  is  applied  to  the  beam  model  on 
element  1. 

2.  The  mass  matrix  is  assembled  for  the  mass  perturbation,  then  the  mass 

matrix  and  partial  derivative  are  calculated. 
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3.  The  first  column  of  the  sensitivity  matrix  is  calculated  using: 


“  A  At  ' 

r  ^ 

=  {®,} 

ADV  _ 

l  ) 

A K  .  AM 

- A  i ' 


A DV  ADV 


(36) 


Note:  This  equation  is  part  of  derivation  of  the  sensitivity  matrix  that  is 
presented  in  the  next  section. 

4.  The  process  is  now  repeated  for  element  2,  and  continues  to  be 
completed  for  all  elements. 

5.  Once  the  mass  portion  has  been  calculated  the  same  process  is  repeated 
with  respect  to  the  stiffness  of  each  element.  Again  a  1%  perturbation 
was  used. 


6.  The  two  separate  portions  (mass  and  stiffness)  are  then  combined  to 
form  the  total  sensitivity  matrix  shown  below.  For  this  study  only  mass 
and  stiffness  varied,  however  many  other  design  variables  can  be  used  in 
the  sensitivity  matrix. 


PI 


dwi  dwf 
dDv  i  dDv  2 

dwl  dw2n 
dDv  i  dDv  2 


(37) 


This  is  considered  the  base  system  sensitivity  matrix  because  it  does  not  take  into 
account  any  artificial  boundary  conditions.  When  adding  ABC  the  process  is  repeated, 
but  the  mode  shapes  used  are  those  of  the  ABC  system.  The  total  system  is  described 
below  including  n  number  of  ABC. 


22 


dw f 

dwf 

dDv  i 

dDvi 

"dw2 

dDv  i 

dDv  2 

'  dwf 

dwf 

dDv  i 

dDv  2 

dw„2 

dDv  i 

.. 

dDv  2 

dwf 

dwf 

dDv  i 

dDv  2 

dDv  i 

.. 

dDvi 

(38) 


This  shows  that  a  change  in  model  parameters  at  one  element  will  alter  the  modal 
parameters  of  the  entire  model.  However,  the  amount  of  change  in  modal  parameters  is 
dictated  by  the  sensitivity  matrix,  allowing  a  perturbation  to  have  significant  or  minimal 
effect.  It  is  because  of  the  influence  that  the  sensitivity  matrix  is  a  key  component  in 
determining  error  prediction  in  a  finite  element  model.  Several  properties  of  the 
sensitivity  matrix  will  be  investigated  later  to  understand  its’  ability  in  error  detection  and 
prediction 

B.  DERIVATION  OF  SENSITIVITY  MATRIX 

To  generate  the  sensitivity  matrix,  starting  with  the  eigenvalue  problem  of  a 
conservative  n  DOF 

[X-/LM]  {$,}  =  {()}  (39) 

Now  taking  Eqn.  (39)  and  taking  the  derivative  with  respect  to  the  design  variable 
generates  Eqn.  (40). 
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A K 
A DV 


.  AM 

A/ - 

A DV 


Ah  ' 
A DV 


{O  i}  +  [K-AM\ 


AO/  1 
A DV* 


=  {0} 


Expanding  Eqn.  (40)  and  pre-multiplying  by  {O}  leads  to  Eqn.  (41) 


(40) 


O 


"  A K  " 

/  ■\  *  /  \T 

“  AM" 

I®!- 

“  Ah  " 

ADV 

{Oj-/L{Oj 

ADV 

ADV  _ 

Now  invoking  the  use  of  the  matrix  identity  {a}7  [/>]  {c}  =  { c}7  [/)]  [<;/ J  the  last  term 

on  the  left  hand  side  can  be  replaced  by  Eqn.  (42),  which  applying  Eqn.  (39)  will  lead  to 
Eqn.  (43) 


(42) 


'  AO/  j7 

adv\ 


[if  -/tM]{0/}  =  0 


(43) 


This  then  reduces  Eqn.  (41)  to 


O; 


A K 
ADV 


{0/}-A«{0/}7 


AM 

ADV 


{O/}  - 


Alt 

ADV 


{0/}r[M]{0/}  =  {0}  (44) 


Now  using  orthogonality,  where  [O/]7  [/V/][0,]  =  1  Eqn.  (44)  further  reduces  to 


O, 


A K 
ADV 


{0/}-/L/{0/}7 


AM 

ADV 


{©,}- 


Ah 

ADV 


(45) 


Now  bringing  the  mass  and  stiffness  portion  to  the  right  hand  side  hand  side 
Ah 


ADV 


=  {0'j 


Aif  AM 
ADV  '  ADV 


{O,} 


(46) 


This  can  be  broken  down  further  into  its  separate  parts  consisting  of  the  mass 
portion  and  stiffness  portion. 
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“  A  At  " 
A DV  _ 

“  A/1,-  “ 
ADV  _ 

A K 
A DV 

.  AM 
— Asi - 


A DV 


{O/}  where  [A A-]  =  [/A- 

{$,}  where  [AM]  =  [Mv-M/] 


(47) 

(48) 


For  both  the  mass  and  stiffness  the  influence  of  the  mode  shapes  is  evident. 
Normalizing  each  row  of  the  mass  and  stiffness  sensitivity  matrix  demonstrates  this. 
Figures  (10)  and  (11)  show  the  base  system  mass  and  stiffness  sensitivity  matrix  rows 
nonnalized  for  the  first  five  rows.  These  first  five  rows  correspond  to  the  first  five  modes 
of  the  base  system.  Each  element  is  represented  by  a  bar  indicating  its  influence  on  the 
change  of  the  systems  natural  frequency  with  a  respective  change  of  design  variable.  For 
instance,  the  stiffness  sensitivity  matrix  indicates  that  a  change  in  stiffness  closer  to  the 
clamped  end  of  the  beam  would  have  a  much  greater  influence  on  the  change  of  natural 
frequency  of  the  system  as  opposed  to  a  change  of  stiffness  at  the  free  end  of  the  beam. 
Just  the  opposite  is  true  for  the  mass  sensitivity  matrix.  It  should  be  noted  that  the  mass 
sensitivity  rows  are  an  absolute  value  and  that  an  increase  of  mass  would  generally  lower 
a  systems  naturally  frequency,  however,  it  was  desired  to  show  the  influence  with  respect 
to  element  position  and  thus  all  magnitudes  were  made  positive. 
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Figure  10.  Base  system,  stiffness  sensitivity  matrix,  modes  1:5 
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Figure  11.  Base  system,  mass  sensitivity  matrix,  modes  1:5 
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From  these  figures  it  can  be  seen  that  where  the  graphs  show  elements  of  low 
sensitivity  correspond  to  areas  where  the  sensitivity  matrix  would  not  identify  a  change  in 
design  variable.  It  should  be  noted  that  due  to  the  normalization  of  each  row  the  graphs 
are  misleading  in  the  frequency  influence.  The  higher  modes  have  a  much  larger  change 
of  frequency  with  a  change  in  design  variable.  In  other  words  a  change  in  the  design 
variable  of  row  one  at  the  element  with  magnitude  of  one  would  produce  a  much  smaller 
change  in  frequency  than  a  change  in  the  design  variable  at  the  element  with  magnitude 
of  one  in  row  two  and  so  on. 

Figures  (12)  and  (13)  show  how  the  application  of  an  ABC  alters  the  respective 
rows  the  sensitivity  matrix.  The  green  triangle  indicates  the  DOF  where  the  ABC  has 
been  applied  to  the  system.  In  general  for  stiffness  the  application  of  an  ABC  drives  up 
the  values  of  the  sensitivity  matrix  in  that  region  where  as  for  mass  it  pushes  the  larger 
values  away  from  the  restrained  DOF.  Again  this  was  intuitive  from  the  resulting  mode 
shapes  when  the  ABC  are  applied.  Figure  (14)  plots  the  modes  for  both  the  base  and 
ABC  system.  The  effect  of  the  ABC  on  the  mode  shapes  corresponds  to  that  of  the 
sensitivity  matrix. 
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Figure  12.  ABC  system,  DOF  31  pinned,  mass  sensitivity  matrix,  modes  1:5 
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Figure  13.  ABC  system,  DOF  31  pinned,  stiffness  sensitivity  matrix,  modes  1:5 
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MODE  SHAPES  FOR  BASE  SYSTEM,  MODES  1:5 


MODE  SHAPES  FOR  ABC  SYSTEM,  MODES  1:5 


Figure  14.  Plot  A.  Base  System,  Modes  1:5.  Plot  B.  ABC  System  DOF  31  pinned  , 

Modes  1:5 

The  relationship  between  mode  shapes  and  the  sensitivity  matrix  often  plays  a  key 
role  in  the  updating  of  a  FEM.  The  use  of  ABC  can  play  a  key  role  in  helping  to  produce 
a  well  conditioned  problem  and  solution  in  the  updating  of  the  model.  As  discussed 
before  a  poorly  conditioned  problem  will  be  one  that  has  a  poor  solution.  However,  a 
well  conditioned  problem,  while  it  is  not  guaranteed  to  have  a  quality  solution  it  is  not 
ruled  out  by  conditioning. 
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VI.  ERROR  PREDICTION 


An  analysis  of  error  prediction  was  conducted  using  MATLAB  to  create  two 
separate  beams.  Each  beam  was  identical  in  dimensions,  42  inches,  1.5  inches,  linch  and 
equal  density  and  Young’s  Modulus.  Each  beam  was  modeled  with  a  twenty  element 
F.E.M.  However,  one  beam  was  considered  the  analytical  beam,  which  represented  the 
baseline  beam,  or  the  F.E.M.  The  second  beam  was  considered  the  experimental  beam 
that  had  a  known  perturbation  applied  at  a  known  element.  Although,  this  beam  too  was 
a  F.E.M  it  represented  an  experimental  structure  from  the  laboratory.  The  errors  imposed 
on  the  experimental  beam  consisted  of  mass  or  stiffness  errors,  generally  in  the  amount  of 
10  percent.  A  diagram  of  the  cantilever  beam  is  shown  below  in  Figure  (15). 


Figure  15.  Cantilever  Beam 


Numerous  trials  were  run  with  the  program  to  help  identify  error  prediction  and 
any  correlations  or  existing  indicators  of  a  strong  or  weak  solution.  Then  for  the  base  and 
all  ABC  error  prediction  was  computed.  This  was  done  for  summations  of  modes  one 
through  twenty.  For  each  of  these  runs  the  condition  of  the  sensitivity  matrix  was 
computed  as  well  as  the  rank.  One  of  the  areas  of  interested  was  the  number  of  modes 
required  to  attain  a  solution  within  ten  percent  of  the  actual  error  at  the  affected  element. 

To  solve  for  the  error  the  matrix  equation  j  Aw2  j  =  [r]  {  ADv}  was  used.  In  order 
to  solve  for  the  change  in  design  variable  { ADv}  MATFAB  arranged  the  equation  as 
follows 

{ADv}  =  \T  ]  1 1  Aw2}  (49) 
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this  is  a  challenging  problem  due  to  the  nature  of  [r]  which  is  not  symmetric  when 

dealing  with  less  than  20  modes  in  the  solution.  Due  to  this  several  different  solution 
methods  were  investigated  to  try  and  develop  better  performance  from  MATLAB.  In 
order  to  fully  understand  the  results  from  the  MATLAB  program  a  brief  discussion  of  the 
solution  methods  can  give  insight  as  to  why  certain  values  were  obtained  or  not  obtained. 

A.  UNDERDETERMINED  MATRIX  ALGEBRA  SOLUTIONS 

Looking  at  the  matrix  problem  is  standard  form 

[A]{x}  =  {b}  (50) 

where  [  A ]  represents  the  sensitivity  matrix,  j  b  j  is  the  square  of  the  change  in  natural 
frequency,  and  jx}  is  the  “error”  or  difference  in  stiffness  or  mass  between  the 
experimental  and  analytical  model.  Or  more  completely 


Where  n  columns  represents  the  elements  of  the  model  and  m  rows  represent  the 
number  of  modes  retained  in  the  sensitivity  matrix.  The  {x}  vector  is  always  of  n  length 

and  the  {b}  vector  is  of  length  m.  With  the  particular  problem  being  addressed  in  <  n 
thus  there  are  several  methods  to  solve  the  problem. 

1.  Condition  and  Rank  of  a  Matrix 

Of  great  importance  in  the  matrix  solution  of  a  problem  is  the  “health”  of  the 
matrices  involved.  Some  of  the  most  important  indicators  of  the  health  of  a  matrix  are  its 
rank  and  condition.  Matrix  rank  can  be  thought  of  as  the  common  dimension  of  the  row 
space  and  column  space  of  that  given  matrix  (Anton,  2005).  Essentially  this  is  saying 
the  rank  is  the  number  of  linearly  independent  rows  that  exist  in  a  matrix.  Thus  if  a  n  by 
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m  matrix  has  a  rank  less  than  n  it  is  considered  singular  and  all  of  its  rows  are  not  linearly 
independent  resulting  in  fewer  solution  equations. 

The  condition  number  of  a  matrix  is  a  measure  of  the  sensitivity  of  the  linear 
system [A]  {x}  =  {£}  .  Essentially  this  is  a  measure  of  the  change  in  the  solution  (“x”) 

with  a  small  change  in  either  “A”  or  “b”  or  both.  A  healthy  matrix  would  have  a 
condition  number  close  to  one,  where  as  a  very  large  condition  number  would  indicate  a 
matrix  that  would  not  lead  to  a  reliable  solution  (Watkins,  2002).  A  good  condition 
number  can  lead  to  a  reliable  solution  depending  on  the  solution  method  that  was  utilized. 
Several  solution  methods  are  discussed  next. 

2.  Pseudoinverse  Method 

The  inverse  of  [^4]  exists  only  if  [d]  is  square  and  has  full  rank.  If  this  is  the 
case  the  solution  is 

w=[^r>}  (52) 

The  pseudoinverse,  [d]  is  a  generalization  of  the  inverse  and  exists  for  any 
matrix  and  the  solution  to  [A]  {x}  =  [b\  becomes 

{x}  =  [zl]>}  (53) 

The  best  way  to  generate  [A]+  is  by  use  of  singular  value  decomposition  which 

provides  a  numerically  robust  solution  to  the  least  squares  problem. 

A  =  USVT  (54) 

Where  U  and  V  are  orthogonal  and  are  ( n,n )  and  S  is  diagonal  of  dimensions 
(, m,n )  with  real,  non-negative  singular  values.  This  leads  to 

A+  =V(STSy'STUr  (55) 
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However,  in  the  under-determined  case  because  the  rank,  r,  of  [ A ]  is  less  than  n, 

(A7  *  A^j  does  not  exist,  thus  the  program  uses  only  the  first  r  singular  values  reducing 

S  to  an  (r,r)  matrix  and  shrinking  U  and  V  accordingly  (Bock,  1998).  Combining  Eqn. 
(54)  and  (55)  leads  to 

x  =  vs~lUTb  (56) 

When  this  solution  method  was  incorporated  in  the  error  solution  problem  it  gave 
results  that  tended  to  distribute  the  error  location  across  the  beam  rather  than  look  for  a 
point  mass  or  stiffness  solution.  As  a  result  the  error  predicted  at  the  affected  element 
was  well  below  the  actual  value.  This  is  shown  in  Figure  (16)  below.  The  “actual”  error 
is  .1  located  at  element  15,  the  Pseudoinverse  method  distributes  the  error  across  all 
elements.  Thus,  for  this  particular  application,  this  method  was  not  seen  as  a  reliable 
solution. 


Pseudoinverse  Method 


Figure  16.  Pseudoinverse  Method  vs.  QR  Decomposition  Method 
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3.  Orthogonal-triangular  (QR)  Decomposition 


Orthogonal-triangular  decomposition  or  QR  decomposition  expresses  the  matrix 
as  the  product  of  a  real  complex  unitary  matrix,  Q  and  an  upper  triangular  matrix,  R.  The 
Householder  method  of  this  type  allows  for  column  pivoting.  Thus,  as  the  program 
searches  for  a  solution  it  will  select  columns  of  [^4]  that  are  more  orthogonal  to  the  other 
used  columns  and  push  less  orthogonal  columns  further  down  in  the  equation.  This  does 
produce  fairly  accurate  results  with  the  number  of  non-zero  scalars  in  the  {x}  vector 

equal  to  the  rank  of[^4] .  This  gives  results  that  do  not  allow  for  the  examination  of  some 

of  the  properties  of  \A\  due  to  the  change  in  columns.  By  using  a  QR  method  that  does 

not  allow  for  column  pivoting,  the  properties  of  [ A ]  are  retained  and  can  be  examined. 

In  particular  the  condition  of  the  sensitivity  matrix  was  examined.  The  use  of  column 
pivoting  did  improve  the  results  in  tenns  of  number  of  modes  needed  to  locate  the  error. 
Table  (1)  and  (2)  below  show  this.  Table  (1)  shows  the  results  with  the  use  of  column 
pivoting  and  indicates  that  the  error  is  found  with  6  modes,  where  as  Table  (2)  shows  the 
results  without  employing  column  pivoting  and  it  takes  9  modes  to  find  the  error. 
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Modes 

1 

1:2 

1:3 

1:4 

1:5 

1:6 

1:7 

1:8 

1:9 

1:10 

Element 

Perceived  Error  (percent) 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3 

0 

0 

0 

0 

0 

0 

0.0126 

0 

0.006 

0.0186 

4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

5 

0 

0 

0.0352 

0.0155 

0.0221 

-0.001 

0.0066 

-0.004 

0 

0.012 

6 

0 

0 

0 

0 

0 

0 

0 

0.01 

0.0059 

-0.004 

7 

0 

-0.008 

0 

0 

-0.012 

0 

0 

0 

0 

-0.002 

8 

0 

0 

-0.022 

-0.008 

0 

0 

-0.002 

-0.005 

-0.004 

0 

9 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

11 

0 

0 

0 

0 

0.0054 

0.0002 

-0.011 

-0.002 

-0.004 

-0.008 

12 

0 

0 

0 

0 

0 

0 

0 

0.0044 

-0.001 

-0.012 

13 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

14 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

15 

0 

0 

0 

0 

0 

0.1015 

0.0924 

0.0946 

0.0965 

0.102 

16 

0 

0 

0 

0.0387 

0 

0 

0 

0 

-0.003 

-0.013 

17 

0 

0 

0 

0 

0.0374 

-0.001 

-0.001 

0.0043 

0.0035 

0.004 

18 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

19 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

20 

0.0034 

0.0035 

0.0039 

-0.006 

-0.005 

0.0003 

0.001 

-0.008 

-0.000 

0.0009 

Table  1 .  QR  Decomposition  with  column  pivoting 
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Modes 

1 

1:2 

1:3 

1:4 

1:5 

1:6 

1:7 

1:8 

1:9 

1:10 

Element 

Perceived  Error  (percent) 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

0 

0 

-0.001 

0.0001 

3 

0 

0 

0 

0 

0 

0.0478 

0.0443 

0.0599 

0 

0.0187 

4 

0 

0 

0 

0 

0.0228 

0 

-0.021 

-0.037 

-0.001 

-0.053 

5 

0 

0 

0 

0.0425 

0 

0.0031 

0 

0 

0.0005 

0 

6 

0 

0 

0.0165 

0 

0.026 

0 

0.0558 

0.0423 

0 

0.0942 

7 

0 

0 

0 

0 

0 

0.0388 

0 

0 

0 

0 

8 

0 

0 

0 

0 

0 

0 

-0.053 

0 

0 

0 

9 

0 

0 

0 

0 

0 

0 

-0.005 

0 

0 

0 

10 

0 

0 

-0.069 

0 

-0.062 

-0.068 

-0.034 

-0.035 

0 

-0.012 

11 

0 

0 

0 

-0.047 

0 

0 

0 

0 

0 

0 

12 

0 

0 

0 

0 

0 

0 

0 

-0.019 

-0.000 

-0.012 

13 

0 

0 

0 

0 

-0.081 

-0.12 

0 

-0.050 

-0.001 

-0.011 

14 

0 

0 

0 

0 

0 

0 

0 

0 

0.0014 

-0.007 

15 

0 

0 

0 

0 

0 

0 

0 

0 

0.0972 

0 

16 

0 

0 

0 

0 

0 

0 

0 

0.0267 

0 

0.0075 

17 

0 

0.0603 

0 

0.0328 

0 

0 

0 

0 

0.0008 

0.0737 

18 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

19 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

20 

0.0415 

0.0042 

0.0485 

0.0278 

0.0687 

0.0809 

0.0471 

0.048 

0 

0 

Table  2.  QR  Decomposition  without  column  pivoting 


The  nature  of  this  method  proved  to  give  results  that  reflected  the  point  mass  or 
stiffness  at  a  particular  element  rather  than  distributing  the  mass  or  stiffness  across  the 
entire  length  of  the  beam  like  the  Pseudoinverse  method  did.  Once  this  method  found  the 
error  at  the  correct  location  it  was  usually  within  ten  percent  of  the  actual  error.  Of 
concern  with  this  method  was  that  it  would  frequently  pick  up  the  error  in  the  correct 
location  and  amount,  retain  the  solution  for  an  additional  two  or  three  modes,  then  lose 
the  solution  completely  at  the  location  as  an  additional  few  modes  were  added.  Then  the 
solution  was  regained  with  the  addition  of  more  modes.  The  reason  for  this  “losing”  of 
the  solution  lies  in  how  the  program  solves  the  underdetennined  problem  and  looks  for 
the  solution  that  minimizes  the  length  of  the  vector  [^4]  {x}  —  {Z?}  .  Thus,  as  the  program 

was  searching  for  the  “shortest”  solution  to  the  problem  at  times  that  solution  would  not 
include  the  desired  element  error.  This  can  be  seen  in  Table  (3)  where  the  error  in 
element  15  is  identified  after  4  modes,  but  then  lost  as  modes  5  and  6  are  incorporated, 
then  found  again  with  mode  7  added.  Then  the  error  was  lost  for  mode  8  then  identified 
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again  as  mode  9  was  incorporated.  It  was  attempted  to  find  a  predictor  that  would 
indicate  when  the  solution  would  drop  out  the  correct  location.  It  was  thought  that  a 
dramatic  increase  in  condition  or  a  decrease  in  rank  might  indicate  when  the  error  would 
be  dropped,  but  these  did  not  turn  out  to  be  indicative  of  when  the  error  would  be 
dropped.  The  change  in  condition  can  be  shown  from  Table  (3)  where  from  modes  4 
though  9  the  error  is  found  then  lost  several  times.  Throughout  the  process  the  condition 
number  increases  but  not  more  dramatically  when  the  error  is  lost  or  regained. 

Tables  (3)  and  (4)  compare  the  results  of  the  Orthogonal-Triangular 
Decomposition  and  the  Pseudoinverse  method.  The  example  is  a  stiffness  error  at 
element  15,  with  an  ABC  applied  at  DOF  9.  Element  15  is  highlighted  in  red  and 
elements  4  and  5  in  blue  to  indicate  that  their  shared  DOF  is  the  DOF  with  the  applied 
ABC.  Modes  1  through  10  were  summed  and  the  condition  of  the  respective  sensitivity 
matrix  was  taken  for  each  summation.  The  table  clearly  shows  how  the  Orthogonal- 
Triangular  Decomposition  method  finds  the  error  at  the  sum  of  modes  1  through  4,  then 
loses  it  at  as  modes  5  and  6  are  applied.  The  table  also  shows  how  the  Pseudoinverse 
method  spreads  out  the  error  and  does  not  perceive  a  point  error. 
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Modes 

1 

1:2 

1:3 

1:4 

1:5 

1:6 

1:7 

1:8 

1:9 

1:10 

Condition 

1 

42.740 

347.26 

1361.1 

3707.6 

7673.3 

13356 

20088 

35573 

58975 

Element 

Perceived  Error  (percent) 

1 

0 

0 

0 

0 

0 

0 

0.001 

0.0501 

0.0002 

0.0006 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3 

0 

0 

0 

0 

0 

0.087 

0 

0 

0 

0 

4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

5 

0.003 

0.0032 

0.0023 

0.0001 

0.0015 

0.0176 

0.0001 

0 

0 

0 

6 

0 

0 

0 

0 

0 

0 

0 

0.0527 

0.0002 

0 

7 

0 

0 

0 

0 

0 

0 

0 

0.0694 

0 

0.0019 

8 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0019 

9 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0005 

10 

0 

0 

0 

0 

0 

0 

0.0007 

0 

0.0007 

0 

11 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

12 

0 

0 

0.0261 

0.0007 

0.0223 

0.0218 

0.0007 

0.0263 

0.0002 

0.0001 

13 

0 

0.0626 

0 

0 

0 

0 

0 

0 

0 

0 

14 

0 

0 

0 

0 

0.0682 

0.0622 

0 

0.047 

0.0012 

0.0022 

15 

0 

0 

0 

0.0922 

0 

0 

0.0926 

0 

0.0928 

0.0946 

16 

0 

0 

0.0878 

0 

0.0555 

0.056 

0 

0.0586 

0.0001 

0.0007 

17 

0 

0 

0 

0.0006 

0 

0 

0.0005 

0 

0 

0 

18 

0 

0 

0 

0 

0.0222 

0.0378 

0 

0.0704 

0.0001 

0.0008 

19 

0 

0 

0 

0 

0 

0 

0 

0.1375 

0.0008 

0.002 

20 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  3. 

Orthogona 

-triangu 

ar  decomposition 

Modes 

1 

1:2 

1:3 

1:4 

1:5 

1:6 

1:7 

1:8 

1:9 

1:10 

Condition 

1 

42.740 

347.26 

1361.1 

3707.6 

7673.3 

13356 

20088 

35573 

58975 

Element 

Perceived  Error  (percent) 

1 

0.0004 

0.0073 

0.0166 

0.0088 

0.0005 

0.0053 

0.0082 

0.0029 

0.0014 

0.0053 

2 

0.0003 

0.0039 

0.0047 

0.0009 

0 

0.001 

0.0038 

0.0024 

0.0017 

0.0076 

3 

0.0003 

0.0016 

0.0003 

0.0013 

0.0003 

0.0052 

0.0097 

0.0029 

0.0009 

0.002 

4 

0.0002 

0.0003 

0.0019 

0.0051 

0.0005 

0.0037 

0.0021 

0.0007 

0.0013 

0.0088 

5 

0.0002 

0 

0.0063 

0.0068 

0.0002 

0.0004 

0.0051 

0.0044 

0.0019 

0.0028 

6 

0.0002 

0.0005 

0.01 

0.0045 

0 

0.0037 

0.0113 

0.0018 

0.0005 

0.008 

7 

0.0001 

0.0015 

0.0105 

0.001 

0.0002 

0.0064 

0.0029 

0.0018 

0.0024 

0.0036 

8 

0.0001 

0.0027 

0.0078 

0.0005 

0.0005 

0.0023 

0.0042 

0.0043 

0.0005 

0.007 

9 

0.0001 

0.0038 

0.0035 

0.0038 

0.0004 

0.0009 

0.0114 

0.0007 

0.002 

0.0047 

10 

0.0001 

0.0045 

0.0005 

0.0073 

0.0001 

0.0054 

0.0035 

0.0033 

0.0012 

0.0058 

11 

0 

0.0048 

0.0008 

0.0072 

0.0001 

0.0054 

0.0035 

0.0033 

0.0011 

0.0059 

12 

0 

0.0046 

0.0044 

0.0036 

0.0004 

0.0009 

0.0114 

0.0007 

0.002 

0.0047 

13 

0 

0.004 

0.0096 

0.0005 

0.0005 

0.0023 

0.0042 

0.0043 

0.0005 

0.007 

14 

0 

0.0031 

0.0134 

0.0012 

0.0002 

0.0063 

0.0028 

0.0018 

0.0024 

0.0036 

15 

0 

0.0021 

0.0141 

0.0053 

0 

0.0036 

0.0112 

0.0017 

0.0005 

0.008 

16 

0 

0.0012 

0.0114 

0.0087 

0.0003 

0.0004 

0.0049 

0.0043 

0.0019 

0.0027 

17 

0 

0.0006 

0.0068 

0.008 

0.0006 

0.0045 

0.0025 

0.0007 

0.0013 

0.0086 

18 

0 

0.0002 

0.0028 

0.0043 

0.0005 

0.0073 

0.0124 

0.0036 

0.0011 

0.0023 

19 

0 

0 

0.0006 

0.0011 

0.0002 

0.0033 

0.0085 

0.0043 

0.0027 

0.0109 

20 

0 

0 

0 

0.0001 

0 

0.0002 

0.0007 

0.0004 

0.0004 

0.002 

Table  4.  Pseudo  inverse  Method 
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B. 


ERROR  PREDICTION  SOLUTION 


Numerous  trials  were  run  to  help  identify  error  prediction  and  any  correlations  or 
existing  indicators  of  a  strong  or  weak  solution.  For  each  element  a  perturbation  of  mass 
or  stiffness  of  ten  percent  was  applied.  Then  for  the  base  and  all  ABC  error  prediction 
was  computed.  This  was  done  for  summations  of  modes  one  through  twenty.  For  each 
of  these  runs  the  condition  of  the  sensitivity  matrix  was  computed  as  well  as  the  rank. 
The  rank  for  every  run  was  identical  to  the  number  of  modes  retained  and  thus  was  not 
seen  as  a  good  predictor  of  solution  accuracy.  The  condition  number  was  greater  for  the 
systems  that  required  more  modes  to  achieve  the  error  tolerance.  But  again  this  was  not 
seen  as  a  good  predictor  because  the  condition  of  the  sensitivity  matrix  is  independent  of 
error  location. 

The  following  example  demonstrates  the  process  and  some  of  the  key  components 
of  the  results.  Using  the  twenty  element  beam  model  element  15  was  subjected  to  the  ten 
percent  mass  or  stiffness  perturbation.  The  system  was  modeled  under  the  base  condition 
(no  ABC),  ABC  with  DOF  9  pinned  then  ABC  with  DOF  31  pinned.  These  conditions 
are  displayed  in  Figure  (17). 


Figure  17.  Twenty  element  beam,  ABC  system  with  DOF  9  pinned  (1)  and  ABC 

system  with  DOF  3 1  pinned  (2). 


Figures  (18)  show  the  results  for  the  applied  error.  In  the  graph  the  red  bars 
indicate  the  base  system,  blue  indicates  the  system  with  DOF  9  pinned,  and  green  with 
DOF  3 1  pinned.  The  stem  plot  indicates  the  magnitude  and  where  the  error  or  damage  is 
located. 
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MODE  1 


MODES  1:6 


Figure  18.  Element  15,  Stiffness 

The  stiffness  perturbation  applied  to  element  1 5  was  difficult  for  the  base  system 
to  pick  out  due  to  its  location  in  a  region  of  the  beam  that  was  relatively  insensitive  to  a 
change  in  stiffness.  The  two  ABC  did  have  better  success  in  detecting  the  error.  In 
particular  the  ABC  applied  at  DOF  3 1  picked  up  the  error  quickly  and  retained  it.  This 
demonstrates  one  of  the  key  findings  in  how  when  trying  to  predict  the  error  location  the 
model  was  highly  influenced  by  areas  of  strong  sensitivity.  Errors  in  those  areas  were 
detected  with  relative  ease.  In  particular  with  stiffness  applying  an  ABC  adjacent  to  the 
element  with  error  can  detect  and  keep  the  error  location  very  well.  This  was  a  reflection 
of  the  increased  sensitivity  when  applying  boundary  conditions  which  had  the  effect  of 
increasing  the  sensitivity  in  the  general  region.  This  can  be  seen  in  Figure  (19)  below 
where  the  top  graph  is  the  base  system,  the  second  ABC  with  DOF  9  pinned,  the  third 
ABC  with  DOF  19  pinned  and  fourth  ABC  with  DOF  29  pinned.  The  one  area  where  this 
was  not  particularly  true  was  at  the  free  end  of  the  beam  where  the  application  of  an  ABC 
did  not  dramatically  increase  the  sensitivity.  This  can  be  seen  as  a  result  of  the  strain 
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energy  at  the  end  of  a  pinned  beam  being  relatively  small  thus  a  change  in  stiffness  at  this 
point  would  not  have  an  strong  influence  on  the  natural  frequency  of  the  system. 


Base  System 
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Figure  19.  Movement  of  pinned  DOF  effect  on  sensitivity 


MODES  1:6 


Figure  20.  Element  15,  Mass 
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When  dealing  with  a  change  of  mass  the  system  tended  to  have  similar  behavior 
but  it  was  somewhat  more  difficult  to  anticipate.  It  was  expected  that  the  base  system 
and  the  ABC  system  with  DOF  9  pinned  would  have  success  locating  the  mass  error  at 
element  15.  This  was  only  true  once  8  or  more  modes  were  retained.  Again  the  system 
performed  better  with  areas  of  high  sensitivity  finding  the  damage  in  the  structure.  This 
was  however,  more  difficult  to  anticipate  with  the  mass  perturbation  due  to  its  behavior 
with  respect  to  additions  of  boundary  conditions. 

Next  an  investigation  of  error  prediction  was  done  with  the  combination  of  base 
and  ABC  modes.  Again  element  15  was  used,  with  ABC  at  DOF  9  and  31.  The  first 
three  modes  for  both  the  base  system  and  ABC  system  were  summed  for  each  of  the  three 
systems  mass  and  stiffness  perturbations.  It  should  be  noted  that  only  one  ABC  at  a  time 
was  combined  with  the  base  system  Figures  (21)  and  (25)  display  the  results.  The  stem 
plot  indicates  the  magnitude  and  location  of  the  error.  The  base  system  in  combination 
with  the  ABC  system  with  DOF  9  pinned  is  displayed  in  blue.  The  base  system  in 
combination  with  the  ABC  system  with  DOF  3 1  pinned  is  displayed  in  green. 


When  dealing  with  the  combination  of  two  boundary  conditions  on  a  system  the 
performance  was  based  on  how  the  boundary  conditions  complimented  each  other  on 
fortifying  areas  of  low  sensitivity.  The  smoother  the  sensitivity  of  the  system,  (fewer 
areas  of  extremely  low  sensitivity)  the  stronger  the  results.  This  can  be  thought  of  as 
combining  sensitivities  of  different  boundary  conditions  to  build  up  areas  of  low 
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sensitivity  to  help  find  the  error.  If  the  error  exists  in  an  area  of  low  sensitivity  then  it 
will  be  difficult  to  locate,  thus  by  eradicating  areas  of  low  sensitivity  error  location  is 
improved.  Figures  (22)  through  (24)  show  the  sensitivities  for  the  base  system,  ABC 
with  DOF  9  pinned  and  ABC  with  DOF  3 1  pinned.  The  pinning  of  both  DOF  9  and 
DOF  3 1  creates  sensitivities  in  areas  that  the  base  system  is  relatively  insensitive.  As  a 
result  the  error  is  detection  is  generally  improved. 


Figure  22.  Base  System,  Modes  1-3,  Stiffness 


Figure  23.  ABC  with  DOF  9  pinned,  Modes  1-3,  Stiffness 


44 


MODE  1 
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Figure  24.  ABC  with  DOF  3 1  pinned,  Modes  1-3,  Stiffness 

Although  the  exact  damaged  element  was  not  found  the  general  region  was  and 
the  magnitude  was  as  well.  The  system  with  ABC  at  Node  3 1  found  the  damage  due  to 
the  high  sensitivity  in  the  region  around  element  15  with  this  particular  boundary 
condition. 


Figure  25.  Element  15,  Base  and  ABC  systems,  Modes  1:3,  Mass 


The  solution  for  a  change  of  mass  demonstrates  this  same  principal.  The  ABC 
applied  further  from  the  clamped  end  of  the  beam  helped  to  even  out  the  sensitivities  of 
the  beam  and  thus  had  improved  error  localization.  The  system  with  DOF  9  pinned 
resulted  in  similar  sensitivities  to  that  of  the  base  system  which  did  not  help  to  strengthen 
the  areas  of  weak  sensitivity.  Figures  (26)  through  (28)  show  the  base  system,  ABC  with 
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DOF  9  pinned  and  ABC  with  DOF  31  pinned.  Of  particular  interest  is  Figure  (28)  and 
how  pinning  DOF  31  lowered  the  sensitivity  around  element  15  and  thus  worsened  the 
error  prediction. 


Figure  26.  Base  System,  Modes  1-3,  Mass 


Figure  27.  ABC  with  DOF  9  pinned,  Modes  1-3,  Mass 
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Figure  28.  ABC  with  DOF  3 1  pinned,  Modes  1-3,  Mass 


For  all  of  the  scenarios  a  pattern  did  emerge  that  depended  on  the  sensitivity 
matrix.  The  magnitude  of  the  perturbation  did  not  have  an  effect  on  error  location  or 
detection.  Rather  the  detection  was  a  function  of  how  sensitive  the  element  was  to  an 
error.  An  element  that  was  located  in  a  region  of  high  sensitivity  resulted  in  relatively 
efficient  error  prediction.  An  element  that  was  located  in  a  region  of  low  sensitivity  had 
a  difficult  time  in  error  prediction. 

An  additional  study  was  done  on  the  stiffness  perturbation.  A  general  solution 
program  procedure  was  proposed.  By  running  the  program  with  all  ABC  applied  to  the 
system  one  at  a  time  it  was  found  that  the  error  could  be  located  by  looking  for  the  largest 
change  in  design  variable  across  the  results.  Figure  (29)  below  shows  the  beam  with  all 
the  possible  ABC  pins.  For  example  a  change  in  stiffness  was  applied  to  element  10,  the 
program  was  then  run  for  all  systems,  base  and  ABC.  Using  just  mode  number  one  it  can 
be  shown  that  use  of  additional  boundary  conditions  can  isolate  the  error.  Figure  (29 
below  shows  where  each  of  the  boundary  conditions  were  applied  one  at  a  time. 
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Figure  29.  Application  of  pinned  ABC  starting  with  clamped  end  of  beam. 


Element 

Base 

ABC  3 

ABC  5 

ABC  7 

ABC  9 

ABC  11 

ABC  13 

ABC  15 

ABC  17 

ABC  19 

1 

0.0133 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

0.0158 

0 

0 

0 

0 

0 

0 

0 

0 

3 

0 

0 

0.0189 

0 

0 

0 

0 

0 

0 

0 

4 

0 

0 

0 

0.0228 

0 

0 

0 

0 

0 

0 

5 

0 

0 

0 

0 

0.0278 

0 

0 

0 

0 

0 

6 

0 

0 

0 

0 

0 

0.0344 

0 

0 

0 

0 

7 

0 

0 

0 

0 

0 

0 

0.043 

0 

0 

0 

8 

0 

0 

0 

0 

0 

0 

0 

0.0545 

0 

0 

9 

0 

0 

0 

0 

0 

0 

0 

0 

0.0703 

0 

10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0924 

11 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

12 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

13 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

14 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

15 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

16 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

17 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

18 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

19 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

20 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Element 

ABC  21 

<\BC  23 

ABC  25 

ABC  27 

ABC  29 

ABC  31 

ABC  33 

ABC  35 

ABC  37 

ABC  39 

ABC  41 

1 

0 

0 

0 

0 

0 

0.0654 

0.0543 

0.0486 

0.0437 

0.0387 

0.0336 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

C 

0 

3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

5 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

6 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

7 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

8 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

9 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

11 

0.0923 

0.0689 

0 

0 

0 

0 

0 

0 

0 

c 

0 

12 

0 

0 

0.0547 

0 

0 

0 

0 

0 

0 

c 

0 
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Table  5.  Error  prediction  using  one  mode  of  each  ABC  system. 


As  the  ABC  approaches  the  error  the  value  of  the  change  in  design  variable 
increases,  then  as  the  ABC  passes  the  value  decreases.  With  this  method  the  error  was 
located  and  the  value  predicted  to  within  7  percent  with  the  use  of  only  one  mode.  The 
behavior  of  the  system  at  the  end  of  the  beam  created  the  poor  prediction  of  the  error 
after  ABC  29,  but  did  not  affect  the  solution  for  errors  located  in  other  regions.  In  the 
case  where  the  error  was  located  in  this  region  this  solution  method  would  not  be  an 
affective  tool.  A  solution  method  of  this  type  was  not  found  to  be  applicable  to  mass 
perturbations. 

C.  OPTIMIZATION  PROGRAM 


Due  to  the  relative  difficulty  of  predicting  mass  error  in  the  structure  an 
optimization  program  was  developed  that  utilized  both  frequency  change  and  mode  shape 
comparison  to  help  locate  the  error.  The  mode  shape  data  was  assessed  using  the  modal 
assurance  criterion  (MAC)(Allemang,  1982). 


MAC 


(57) 


where  | |  is  the  modal  amplitude  at  location  1  of  the  zth  experimental  mode;  | O " y  |  the 

modal  amplitude  at  location  1  of  the  yth  calculated  mode;  and  p  the  length  of  the  modal 
shape  vector  (Zhang,  2000).  The  experimental  eigenvalues  and  eigenvectors  were 
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extracted  from  the  FRFs  taken  from  the  experimental  set  up  described  in  Appendix  A 
which  matches  the  FEM  created  for  the  previous  error  prediction  calculations. 

The  frequency  difference  and  mode  shape  data  were  combined  in  the  objective 
function  below. 
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The  MAC  produces  values  from  0  to  1,  1  indicating  the  two  mode  shapes  being 
compared  are  quite  close.  Thus,  the  MAC  was  subtracted  from  one  to  produce  a  result 
that  would  indicate  a  better  performance  when  minimized.  The  function  was  minimized 
by  use  of  the  MATLAB  algorithm  “fmincom”.  This  program  allows  the  use  of  lower  and 
upper  bounds,  inequalities  and  equalities.  Of  interest  was  the  use  of  just  the  mode  shapes 
or  frequencies  of  the  system  and  if  the  combination  of  the  two  would  produce  a  more 
reliable  and  accurate  solution.  In  addition  a  variety  of  constraints  were  used  in 
performing  the  computations  to  try  and  improve  the  performance  of  the  objective 
function. 


The  problem  was  run  with  no  constraints  and  then  constrained.  The  initial 
constraints  used  in  the  problem  were  the  amount  of  mass  that  could  be  added  to  any  one 
element  as  well  as  the  total  amount  of  mass  that  could  be  added  to  the  system.  Further 
constraints  that  were  recognized  were  the  limitation  of  mass  addition  to  only  certain 
elements. 


1.  Procedure 

To  validate  the  optimization  program  initial  trials  were  run  with  comparing  two 
MATLAB  beams,  one  with  an  error,  the  other  with  no  error.  The  error  was  applied  to 
element  35  of  the  model.  The  magnitude  of  the  error  was  10  percent  of  the  total  beam 
weight.  Three  separate  objective  functions  were  used  in  the  comparison.  The  first  was 
just  the  comparison  of  natural  frequencies, 
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the  second  a  comparison  of  mode  shapes  (MAC), 
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The  above  objective  functions  were  minimized  subject  to  the  following 
constraints. 


1 .  Unconstrained 


or 

2.  Limed  amount  of  mass  that  could  be  added  to  any  one  element  (10  percent), 
or 

3.  Limited  total  amount  of  mass  that  could  be  added  to  the  entire  beam  (10 
percent). 

Once  the  MATLAB  comparison  was  done  the  same  program  was  utilized  in 
comparing  experimental  data  from  the  laboratory.  The  procedure  for  extracting  the 
modal  parameters  from  the  beam  follows. 

The  FRFs  were  taken  from  the  laboratory  beam  at  DOF  42.  This  data  was 
considered  the  base  line  data.  Each  DOF  was  one  inch  apart  and  DOF  1  was  at  the 
clamped  end  of  the  beam.  The  natural  frequencies  and  mode  shapes  extracted  using  the 

RT  Pro  Focus  and  ME’Scope’VES  curve  fit  the  data.  The  curve  fitting  procedure  uses  a 
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mathematical  algorithm  to  estimate  the  modal  parameters  from  the  measured  data.  The 
mode  shapes  extracted  from  the  laboratory  software  consisted  of  real  and  imaginary  data. 
An  algorithm  in  the  software  converted  this  data  to  real  mode  shapes  that  can  be 
compared  to  the  FEM  data  produced  in  MATLAB.  The  algorithm  takes  the  sum  of  the 
squares  of  the  real  and  imaginary  portions  for  the  magnitude,  then  assigns  a  phase  angle 
of  either  0  or  180  degrees.  This  is  also  known  as  the  “Simple  Method”  (Ewins,  2005) 

There  was  an  initial  difference  in  natural  frequencies  between  the  laboratory  beam 
and  the  FEM  that  varied  between  approximately  three  and  five  percent.  Due  to  this 
difference  the  optimization  program  was  run  to  update  the  FEM  to  match  the  natural 
frequencies  of  the  actual  beam  set  up. 

Next  a  mass  error  was  added  to  the  laboratory  beam  at  a  particular  element.  The 
FRFs  were  taken  for  all  42  DOF  and  natural  frequencies  and  mode  shapes  extracted.  The 
optimization  program  was  then  run  again  with  the  updated  FEM  to  locate  the  mass  error. 

2.  Results 

Full  results  are  available  in  Appendix  B. 

1.  Natural  Frequency  only:  The  natural  frequency  used  alone  did  not 
provide  highly  accurate  error  predictions.  The  best  error  prediction  it 
could  accomplish  was  within  74  percent  of  the  actual  error.  For 
several  of  the  cases  the  number  of  iterations  were  exceeded  so  the 
program  stopped  prior  to  meeting  the  desired  tolerance. 

2.  Mode  shapes  only:  Overall,  the  best  performance  was  from  the  use  of 
the  MAC  (mode  shapes  alone).  This  was  the  case  for  both 
experimental  and  MATLAB  data.  Additionally,  as  the  problem 
became  more  constrained  the  error  prediction  became  better.  The  best 
performance  for  this  was  the  fully  constrained  problem  with 
constraints  on  both  total  mass  added  to  the  system  and  total  mass 
allowed  for  each  element.  This  gave  error  predictions  within  5  percent 
for  the  MATLAB  portion  and  9  percent  for  the  experimental. 
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3. 


Mode  shapes  and  natural  frequencies:  The  use  of  both  modes  shapes 
and  natural  frequencies  in  the  objective  function  did  improve  the 
results  from  just  the  use  of  natural  frequencies,  but  did  not  improve  the 
results  from  just  using  mode  shapes  alone. 

The  optimization  program  did  generate  reliable  data  for  the  comparison  of 
eigenvectors  (mode  shapes)  however,  did  not  present  itself  as  a  good  error  predictor  when 
using  only  eigenvalues  (natural  frequencies).  This  was  more  than  likely  due  to  the  fact 
that  when  comparing  eigenvectors  there  is  more  information  about  the  system  being 
compared  (displacement  at  each  DOF).  Where  as  when  comparing  eigenvalues  the  only 
comparison  is  of  several  values,  each  of  which  is  described  by  the  entire  system.  The  use 
of  only  natural  frequencies  in  damage  detection  has  been  known  to  have  limitations 
(Salawu,  1997). 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  use  of  Artificial  Boundary  Conditions  in  sensitivity-based  model  updating 
has  proven  to  be  an  effective  method  of  overcoming  the  shortfalls  of  an  underdetermined 
problem.  The  focus  of  this  thesis  was  to  try  and  evaluate  and  obtain  reliable  methods  to 
help  generate  an  improved  error  location  method. 

A.  CONCLUSIONS 

For  the  system  tested  in  this  thesis,  with  respect  to  perturbations  or  errors  in  mass 
and  stiffness,  the  performance  of  the  updating  could  be  greatly  improved  by  judicious 
application  of  ABC.  For  both  mass  and  stiffness  errors  the  best  results  were  obtained 
when  the  sensitivities  of  the  system  were  more  evenly  distributed  across  the  beam. 
When  a  system  had  areas  of  low  sensitivities  it  created  the  possibility  of  poor  error 
detection.  This  could  be  reduced  by  increasing  the  sensitivity  in  these  regions.  This 
could  be  accomplished  by  applying  an  ABC  in  a  position  that  would  increase  the 
sensitivity  in  the  region(s)  of  poor  sensitivity. 

This  was  relatively  easy  to  anticipate  with  the  lower  modes  of  the  system.  As  the 
modes  increased  so  too  did  the  complexity  of  the  mode  shapes  and  sensitivity 
distribution.  The  sensitivity  distribution  with  respect  to  stiffness  was  particularly  reliable 
to  anticipated,  essentially  knowing  that  the  sensitivity  would  be  increased  in  the  region  of 
the  boundary  conditions  (pin).  The  sensitivity  distributions  with  respect  to  mass  were 
more  difficult  to  predict.  They  tended  to  increase  as  the  distance  from  a  boundary 
condition  increased  which  correspond  to  areas  of  large  displacement  as  these  are  also 
areas  of  large  acceleration  and  hence  kinetic  energy. 

A  better  performing  MATLAB  computation  method  was  not  discovered.  The  QR 
Decomposition  method  produced  results  that  often  dropped  out  the  desired  solution.  A 
correction  for  this  discrepancy  was  not  discovered. 
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The  use  of  an  optimization  program  to  develop  a  solution  method  to  a  mass 
perturbation  was  also  incorporated.  This  program  made  us  of  both  natural  frequencies 
and  mode  shape  information.  This  was  done  in  order  to  improve  the  performance  of  the 
error  prediction. 

B.  RECOMMMENDATIONS 

1.  Investigate  the  properties  of  the  QR  decomposition  algorithm  which  do 
not  smear  the  solution.  This  could  lead  to  an  improved  method  of  error  location  and 
retention  with  addition  modes  summed. 

2.  Comparison  of  mass  and  stiffness  sensitivity  matrices  with  strain  and 
kinetic  energy. 

3.  Develop  an  optimization  program  that  would  minimize  the  differences  in 
modal  parameters  between  a  FEM  and  an  actual  structure  with  a  stiffness  error  applied  to 
the  actual  structure. 

4.  Further  investigate  the  use  of  an  optimization  program  to  select  ABC  that 
would  minimize  areas  of  low  sensitivity  in  a  FEM  to  correlate  to  accurate  model 
updating. 

5.  Test  the  use  of  sensitivity-based  updating  on  a  more  complex  base  system 
(additional  base  boundary  conditions)  to  evaluate  the  minimizing  of  regions  of  low 
sensitivity. 
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APPENDIX  A 


A.  EXPERIMENTAL  SET  UP 

A  block  of  steel  1 8  inches  in  length,  8  inches  wide  and  2  inches  thick  was  placed 
on  a  platform  as  a  foundation.  A  cantilever  beam  made  of  T-6061  Aluminum,  48  inches 

3 

in  length,  1.5  inches  wide  and  0.5  inches  thick,  density  of  0.11  lbf/in  and  elasticity 

2 

modulus  of  10  E6  lbf/sec  -in  was  placed  on  top  of  foundation  steel  table.  The  beam 
extended  42  inches  and  was  clamped  to  the  foundation  with  two  large  C  clamps  and  two 
shorter  length  steel  beams,  each  6  inches  long,  2  inches  wide  and  1  inch  thick.  These 
assisted  in  the  elimination  of  platform  generated  modes. 


Figure  30.  Experimental  beam  set  up 

The  beam  consisted  of  42  elements,  each  1  inch  in  length;  this  corresponded  to  FE 
model  element  quantity  and  length.  A  Series  336  FLEXCEL  ICP  accelerometer  (serial 
number  10860)  was  threaded  into  position  at  node  41  of  the  beam  and  wired  into  Channel 
2  on  a  DACTRON  Focus  front  end  digital  signal  processor  (DSP).  An  excitation  was 
applied  by  a  PCB  Series  086  B03  impact  hammer  (serial  number  269),  which  was  wired 
into  Channel  1  on  the  DSP.  The  accelerometer  and  force  hammer  was  calibrated  using  a 
pendulous  test  and  the  sensitivity  adjustment  was  applied  in  the  set-up  of  RT  Pro  Focus 
5.57  software. 
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B.  DATA  COLLECTION 

Using  DACTRON  RT  Pro  Focus  5.57  software,  FRFs  were  collected  when  the 
roving  force  was  applied  by  the  load  cell  at  each  node.  Node  41  remained  the  reference 
as  the  load  cell  roved  from  one  node  to  next  allowing  for  the  measurement  of  the 
response  at  each  node.  Due  to  this  set-up  only  one  column  or  row  of  the  complete  FRF 
matrix  [H]  was  actually  measured. 

1.  A  RT  Pro  Focus  5.57  software  “Real-time”  project  was  configured  to  measure 
3200  spectral  lines,  8192  points,  with  a  delta  T  of  166. 7ps  over  the  frequency  range  0- 
2400Hz.  The  frequency  range  of  0-2400  Hz  was  chosen  because  it  covered  the  first  10 
modes  of  system  and  signal  resolution  was  sufficient  for  data  acquisition.  The  excitation 
signal  proved  to  be  clean  and  thus  no  window  was  used  for  data  measuring. 

2.  Channel  set-up 

Channel  1  (Excitation)  Channel  2(Response) 

Max  Volts  (mV):  0.1  0.3 

Quantity:  Force  Accel. 

EU:  lbf  gn 

mv/EU:  8.67  104.383 

Coupling:  ICP  AC  0.7  Hz  ICP  AC  0.7  Hz 

Sensitivity  Adjustment:  0  0 


3.  Trigger  set-up 

Source:  Analog  input 

Run  Mode:  Manual  Ann  every  frame 

Input:  Channel  1 

Slope:  Bi-polar 
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Level  (%):  1,  Level  (V):  0 
Pre/Post  Points  (-/+):  -10 
Pre/Post  Time  (-/+):  -1 .67jo,s 

4.  Average  set-up: 

Type:  Linear 
Domain:  Frequency 

Frames:  3  (Each  node  was  excited  3  times  and  an  average  taken  and  saved.) 
Accept/Reject:  Manual  Accept/Reject  every  frame. 

(The  user  rejects  double  taps,  under  powered  or  overloaded  signals.) 

5.  Modal  Coordinate  set-up: 

Auto  increment:  ON 
Rove:  Excitation 

Point  increment:  1 

Export:  UFF  text  fonnat,  frequency  response. 
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APPENDIX  B 


Natural  Frequencies 

MATLAB 

Experimental 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Iterations 

67 

78 

24 

66 

67 

66 

Fct.  Ct 

4242 

4217 

1321 

4224 

3490 

4242 

Element 

1 

0.2098 

0.0991 

0 

-0.0984 

-0.1147 

-0.1147 

2 

0 

0.0732 

0.0009 

0.67 

0.6608 

0.8426 

3 

0 

0 

0 

-0.1886 

-0.1886 

-0.1886 

4 

0 

0 

0 

0.1173 

0.1785 

0.1173 

5 

0 

0 

0 

0.0624 

0.0624 

0.0679 

6 

0 

0 

0.0374 

0.0066 

0.0066 

0.0975 

7 

0.278 

0.0695 

0.0882 

-0.0286 

-0.0286 

-0.0286 

8 

0.0932 

0.2259 

0 

0.0933 

0.0724 

0.1614 

9 

0.004 

0.0171 

0 

-0.0405 

-0.047 

-0.0156 

10 

0.15 

0.0427 

0 

0.1819 

0.1819 

0.1819 

11 

0.0169 

0.0447 

0 

0.1513 

0.1513 

0.1513 

12 

0 

0 

0 

0.1409 

0.1409 

0.1409 

13 

0 

0 

0.035 

-0.0506 

-0.0506 

-0.0506 

14 

0 

0 

0 

0.1444 

0.1444 

0.1444 

15 

0 

0 

0 

-0.0208 

-0.0208 

-0.0208 

16 

0.0013 

0 

0 

0.184 

0.0802 

0.0802 

17 

0.2109 

0.035 

0 

0.1242 

0.1242 

0.1242 

18 

0 

0 

0 

0.0274 

0.0274 

0.0274 

19 

0 

0 

0 

0.0584 

0.0584 

0.0584 

20 

0.013 

0 

0 

-0.0352 

-0.0352 

-0.0352 

21 

0.0252 

0 

0.0153 

-0.0103 

-0.0103 

-0.0103 

22 

0 

0 

0 

0.1608 

0.157 

0.1445 

23 

0.0024 

0 

0 

0.0431 

0.0431 

0.0431 

24 

0.1112 

0 

0.048 

0.0821 

0.0821 

0.0907 

25 

0 

0.0897 

0 

0.0663 

0.0663 

0.0663 

26 

0 

0 

0.039 

-0.0768 

-0.0768 

-0.0768 

27 

0.0003 

0.0432 

0 

0.6871 

0.6015 

0.6939 

28 

0 

0 

0.0352 

0.1137 

0.0581 

0.0581 

29 

0.0479 

0.0018 

0 

0.2309 

0.2309 

0.2309 

30 

0.1194 

0.2332 

0 

0.2467 

0.1548 

0.1034 

31 

0.0007 

0.0387 

0 

-0.1399 

-0.1399 

-0.1399 

32 

0.0604 

0.1613 

0 

-0.0449 

-0.0449 

0.0274 

33 

0.1188 

0 

0 

0.2998 

0.3038 

0.2633 

34 

0.7715 

0.9104 

0.3767 

0.1642 

0.144 

0.2326 

35 

1.4782 

1.539 

0.2024 

0.2385 

0.3174 

0.0821 

36 

0.0604 

0.0793 

0.3749 

0.6164 

0.7565 

0.6955 

37 

0 

0.0212 

0.213 

0.0312 

0.0312 

0.0312 

38 

0.0663 

0 

0.1581 

0.0444 

0.0444 

0.0444 

39 

0.6677 

0.7061 

0 

-0.0786 

-0.0011 

-0.0782 

40 

0.0459 

0.0024 

0 

0.0186 

0.1087 

0.0186 

41 

0.0008 

0.1152 

0.4847 

0.084 

0.084 

0.084 

42 

0 

0 

0.2686 

0.2984 

0.1897 

0.2395 

Table  6.  Objective  Function  with  Natural  Frequencies 
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Mode  Shapes 

MATLAB 

Experimental 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Element 

1 

0.0046 

0.0046 

0.0002 

0.1565 

0.3236 

0 

2 

0.0416 

0.0416 

0 

0.571 

0.5263 

0 

3 

0.1828 

0.1828 

IE-04 

1.2007 

1.0048 

0 

4 

0.3305 

0.3305 

0.0778 

1.6021 

1.1174 

0 

5 

0.3206 

0.3206 

0.0028 

0.3853 

0.3414 

0.2939 

6 

0.1107 

0.1107 

0.0162 

0.4479 

0.0458 

0 

7 

0.2706 

0.2706 

0.0144 

1.0487 

0.0882 

0 

8 

0.3211 

0.3211 

0.0237 

1.8214 

0.7641 

0 

9 

0.2001 

0.2001 

0.0164 

1.6041 

1.0079 

0.2646 

10 

0.2068 

0.2068 

0.0685 

1.1049 

0.867 

0 

11 

0.315 

0.315 

0 

0.1776 

0.0057 

0 

12 

0.2615 

0.2615 

0.0105 

1.3055 

0.3166 

0 

13 

0.2157 

0.2157 

0.0123 

0.2975 

0.0917 

0 

14 

0.2587 

0.2587 

0.0028 

0.2282 

0.0032 

0 

15 

0.2879 

0.2879 

0.0269 

0.83 

0.4894 

0 

16 

0.2654 

0.2654 

0.0752 

1.0372 

0.4838 

0 

17 

0.2133 

0.2133 

0 

0.1547 

0.0229 

0 

18 

0.1869 

0.1869 

0 

0.4684 

0.035 

0 

19 

0.2529 

0.2529 

0.0039 

0.3138 

0.0193 

0 

20 

0.3319 

0.3319 

0.0848 

0.9953 

0.2707 

0 

21 

0.2443 

0.2443 

0 

1.0559 

0.5797 

0 

22 

0.1733 

0.1733 

0.012 

0.3945 

0.3839 

0 

23 

0.2353 

0.2353 

0.0449 

0.317 

0.0012 

0 

24 

0.3023 

0.3023 

0.0348 

0.6827 

0.0313 

0 

25 

0.3175 

0.3175 

0.0063 

0.8176 

0.2061 

0 

26 

0.2266 

0.2266 

0.0032 

0.5868 

0.2703 

0 

27 

0.0536 

0.0536 

0 

0.1869 

0.1067 

0 

28 

0.2986 

0.2986 

0.0859 

0.1289 

0.1753 

0 

29 

0.4143 

0.4143 

0.0123 

0.8184 

0.0357 

0 

30 

0.2596 

0.2596 

0.0021 

0.5693 

0.0059 

0 

31 

0.1672 

0.1672 

0.0022 

0.2602 

0 

0 

32 

0 

0 

0.0004 

0.0071 

0 

0 

33 

0.0986 

0.0986 

0.0007 

0 

0 

0 

34 

2.0037 

2.0037 

0.2836 

4.8772 

3.4007 

0 

35 

2.8639 

2.8639 

3.9957 

10.7002 

4.2 

3.819 

36 

1.8299 

1.8299 

0.1881 

8.6623 

4.2 

1.1225 

37 

0.0407 

0.0407 

0.0043 

0.4437 

0.8355 

0 

38 

0.0482 

0.0482 

0.019 

0.0352 

0 

0 

39 

0.2049 

0.2049 

0.0098 

0.0074 

0 

0 

40 

0.2846 

0.2846 

0.0166 

0.2684 

0 

0 

41 

0.4631 

0.4631 

0.0132 

1.7755 

0.4011 

0 

42 

0.1444 

0.1444 

0.0266 

0.3576 

0.3763 

0 

Table  7.  Objective  Function  with  Mode  Shapes 
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Natural  Frequencies  and  Mode  Shapes 

MATLAB 

Experimental 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Unconstrained 

Constrained 
Element  Mass 

Constrained 
Total  Mass 

Iterations 

78 

79 

34 

66 

72 

39 

Fct.  Ct 

4242 

4217 

1806 

2905 

3247 

1734 

Element 

1 

0.1498 

0 

0.2067 

0.053 

0.1121 

0 

2 

0.3398 

0.5238 

0 

0.0492 

0.7151 

0 

3 

0.0263 

0.0218 

0 

1.0797 

0.9828 

0 

4 

0.1714 

0.0049 

0.2189 

0.446 

0.8307 

0.1022 

5 

0 

0 

0.037 

0.9709 

0.4619 

0.1764 

6 

0 

0 

0 

0.6148 

0.0831 

0 

7 

0.0014 

0.2368 

0 

0.4279 

0.0745 

0 

8 

0.0028 

0 

0 

0.5951 

0.6515 

0 

9 

0 

0.056 

0 

2.2734 

1.0106 

0.3224 

10 

0.2066 

0 

0 

0.6441 

0.8483 

0.0031 

11 

0.0049 

0 

0 

0.0183 

0.0215 

0 

12 

0 

0 

0 

0.6911 

0.2321 

0 

13 

0 

0 

0 

0.3131 

0.0585 

0 

14 

0.0099 

0.0012 

0 

0.019 

0.0412 

0 

15 

0 

0.0295 

0 

1.2532 

0.521 

0 

16 

0 

0.0561 

0 

0.0997 

0.3509 

0 

17 

0.0414 

0.0042 

0 

0.2084 

0.0454 

0 

18 

0.0008 

0.0173 

0 

0.2054 

0.0013 

0 

19 

0.0198 

0.4014 

0 

0.1378 

0.0019 

0 

20 

0.0003 

0 

0 

1.2053 

0.4178 

0 

21 

0.0037 

0.1091 

0 

0.162 

0.3559 

0 

22 

0 

0 

0 

0.0983 

0.4215 

0 

23 

0.0125 

0 

0 

0.8814 

0.0105 

0 

24 

0.0226 

0.0586 

0 

0.39 

0.0272 

0 

25 

0.0342 

0 

0 

0.0491 

0.2321 

0 

26 

0.0053 

0.0048 

0 

0.3196 

0.1301 

0 

27 

0.1928 

0 

0 

0.6125 

0.1351 

0 

28 

0.0214 

0.0772 

0 

0.0172 

0.1504 

0 

29 

0 

0.1595 

0 

0.237 

0.0416 

0 

30 

0.1097 

0.0084 

0 

0.4457 

0 

0 

31 

0 

0.0123 

0 

0.015 

0 

0 

32 

0.0082 

0.0012 

0.1026 

0.0104 

0 

0 

33 

0.6786 

0.7911 

0 

0.1034 

0 

0 

34 

0.0061 

0 

0.6801 

3.2499 

3.313 

0 

35 

1.3859 

1.1767 

1.5291 

9.295 

4.2 

3.5997 

36 

0.1608 

0.07 

0.3855 

7.5685 

4.2 

1.2961 

37 

0.4339 

0.5896 

0 

0.1267 

0.6179 

0 

38 

0.0053 

0.0002 

0 

0.0594 

0 

0 

39 

IE-04 

0 

0 

0.0199 

0 

0 

40 

0.9151 

0 

0 

0.1146 

0 

0 

41 

0 

0.1381 

0 

0.9672 

0.2685 

0 

42 

0 

0 

0.0725 

0.4176 

0.4014 

0 

Table  8. 


Objective  Function  with  Natural  Frequencies  and  Mode  Shapes 
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APPENDIX  C 


A.  ABCRUNTHRU  JRM 


%  ********************  ABCrunTHRU  j  r  m.  m  *********************** 

%  This  propram  calculates  the  condition  number  of  the  following 
%  sensitivity  matrices  used  to  calculate  the  DV  (error  prediction). 
%  1)  Base  system  only  5  modes  (underdeter mi  ned) 

%  2)  ABC  system  10  modes 

%  3 )  Base  system  5  modes  +  5  modes  from  ABC  system 
%  The  last  system  is  calculated  3  times.  Once  for  modes  1-5, 

%  another  for  modes  6-10,  and  again  for  modes  11-15. 

% 

%  This  program  is  called  from  Build2Beams.m. 

% 

%  Written  by  Constance  Ft  S  Fernandez,  Spring  2004 
%  Updated  by  John  R  Mentzer,  Spring  2007 

%l  NPUTS 


%  i  c  n  t  _  o  s  e  t 
%  T  s  e  n  s  _  t  o  t 
%  vect_l  am_tot 

%  OUTPUTS 


%  a  a 
%  dv  a 
%  dv" b 
%  d  v  c 
%  i  n  t  e  r  v  e  I 
%  mode 
%  s  t  a  r  t  mo  d  e 
%  ten 

%  d  v  c  a  I  ABC-  ma  t  r  i  x 
%dv"calABCten  -  matrix 
%dv"cal  BasePlus  -  matrix 
%  cond_ABC  -  mat  r  i  x 
%  cond_ABCten  -  mat  r  i  x 
%  cond_basePlus  -  matrix 

% - I  N I  T I  A L I  ZAT I  ON . 

a  be  =  0; 
ten  =  1 ; 
i  n  t  e  r  v  e  I  =  1; 
i  nt  abc  =  1; 

for"aa  =  1 :  i  c  n  t  oset  +1  %  number  of  conditions  (base  +  ABC) 
a_c  =  1;  %  reinitialize  for  each  ABC  system 

for  mode  =  1:3  %  3  sets  of  modes  per  ABC  system  (10  element  beam) 
%  (  modes  1:  5,  6:  10,  11:  15) 
start  mode  =  abc  +  a  _  c : 

d v  a  =  [  s  t  a  r  t  mo d e :  s  t  a r  t  mode +4] ;  %  modes 

dv"dl=[ten:ten];  %model 

dv"d2  =  [ten+l:ten+l];  %  mode  2 

dv"d3  =  [  t  en+2:  t  en+2]  ;  %  mode  3 

dv"d4  =  [  t  en+3:  t  en+3] ;  %  mode  4 

dv"d5  =  [  t  en+4:  t  en+4] ;  %  mode  5 


dv  c  =  [ten:  ten +9];  %  the  first  10  modes  of  each  ABC  system 
%  [  mo  d  e  s  1-10  only) 

if  dv  a  ==  [  1:  .  25*si  ze(T_sens  tot, 2)];  %  if  sensitivity  matrix 
%has  only  5  rows  than  the" modes  used  are  all  5,  else  modes 
%used  are  first  five  (base)  and  a  set  of  selected  5  modes 
%of  ABC  system  for  a  total  of  10  modes, 
d v _ b  =  [ dv  a] ; 

else 
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d v _ b  =  [1:  .  25*si  ze(T  sens  tot,  2),  dv  a]; 
end 

% —  Base  System  onl  y  —  (underdetermi  ned) 

%save  DV  calculates  of  as  matrix  for  plotting 

dv  cal  ABC( :  ,  i  ntervel  )  =  T  sens_tot(dv  a,  :  )\vect_l  am  tot(dv  a);  hold  on 
%  condition  number  of  Sensitivity  matrix  used  in  DV  cal, 
cond  ABC(  i  ntervel  ,  1)  =  cond(T  sens_t  ot  (  dv  a,:)); 

%  cond  ABC(intervel.l)  = 

c  o  n  d  ( (  T  sens  tot  ( d  v  a ,  :  ) )  *  (  T  sens  tot(dv  a, :))');  %(Ttrans)(T) 

%  ~  " 

%  %mo  d  e  1 

%  dv  cal  _  A  B  C 1  =  T  sens_t  ot  ( dv  d  1 ,  :  )  \  v  e  c  t  lam  t  o  t  ( d  v  dl); 

%  %mode  2 

%  dv  cal  _ A B C 2  =T  sens_t  ot  ( dv  d 2 ,  :  )  \  v e c t  lam  t  o t  ( d v  d 2 ) ; 

%  %mode  3 

%  dv  cal  _ A B C 3  =T  sens_t  ot  (  dv  d 3 ,  :  )  \  v e c t  lam  t  o t  ( d v  d 3 )  ; 

%  %mode  4“ 

%  dv  cal  ABC4  =  T  sens_t  ot  ( dv  d  4 , : )  \  v  e  c  t  lam  totfdv  d  4 ) ; 

%  %mode  5 

%  dv_cal_ABC5  =  T  sens_t  ot  ( dv  d  5 ,  :  )  \  v  e  c  t  lam  t  o  t ( d  v  d  5 ) ; 


%  —  Base  +  5  modes  of  ABC  system - 

dv  cal  _BasePI  us( : ,  i  ntervel  )  =  T  sens  tot(dv  b ,  :  )  \  v  e  c  t  lam_tot(dv_b); 
cond  basePI  us(i  ntervel  ,  1)  =  cond(T  sens_tot[dv  b ,  :  ) ) ; " 

%  cond  basePI  us(  i  ntervel ,  1)  = 

cond(((T_sens_tot(dv_b,:))*(T_sens_tot(dv_b,:))'));%(Ttrans)(T) 

intervel  =  intervel  +  1; 
a _ c  =  a _ c  +  5 ;  %  f  I  v e  modes  used  at  a  t i  me 

end  %  "  mo  d  e "  loop 

%  ---ABC  system  only  ---- 


dv  cal  ABCtenf :  ,  i  nt_abc)  =  T  sens_t  ot  ( dv  c,  :  )\vect_l  am  t  ot  ( dv  c); 
cond  ABCtenjint  abc,  1)  =  condfT  sens  tot[dv  c ,  :  ) ) ; 

%c o n d  ABCt  en(  i  nt  abc,  1)  = 

cond(((T_sens“tot(dv_c,  :T)*(T_sens_tot(dv_c,  : ))' ));  %(Ttrans)(T) 
i  n t  _ a b c  =  i  nt  abc  +1; 

abc  =  abc  +39;  “/(.advances  to  next  ABC  system,  must  change  number  "19" 
%to  reflect  the  number  of  DOF  in  beam.  This  beam  had  10  elements  thus 
%  19  DOF. 

ten  =  ten  +  39;  %  advances  to  next  ABC  system 
end  %  " a  a "  loop 


o/0  ********************  end  ABCrunTHRU  jrm.  m  *********************** 
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B.  ADDLUMPMASS  JRM 


%  ******** ************  AddLumpMass  jrm.  m  *********************** 

%  This  script  constructs  a  vector  of  I  u mp e d  ma s s e s 
%  which  is  added  to  the  diagonal  of  the  BeamX  mass  matrix. 

%  Mass  added  to  [mx]  in  As  s  emb  I  e2  Bea  ms ,  m 

% 

%  Written  by  Prof  J.H.  Gordis 

%  I  n  p  u  t  s  needed 

% . 

%  n  u  m  e  I  e  me  n  t  s 
%  mx 

%  Out  put  s 

% . 

%  ma  s  s  _  d  i  a  g 
%  ma  s  s  _  n  o  d  e 
%  ma  s  s  _  c  o  o  r  d 
%  mas s_ DOF 
%  mx  ■  updated 

d i  s  p ( '  1  ) ;  di  sp( 1  1  ) ; 

di  sp( 1  *********************************************************  j 
d  i  s  p  ( 1  *  *  *  *  Lumped  mass  addition  to  bea  ms  ****') 

dispj'  ****  Lumpmass  DOFs  defined  for  UNRESTRAINED  beams  ****') 
di  sp( 1  *********************************************************  j 

%  i  n  i  t  a  I  i  z  e 

if  exist('mass_diag')  ==  0;  %  define  and  apply  lumped  mass  vector, 

add  mass  =  1  n 1  ; 

add" mass  =  input)'  Add  lumped  masses  to  BeamX  ?  (y/n)  1  ,'s1); 

%  Initialize  vector  to  add  to  [mx]  diagonal, 
ma s s _ d i  a g  =  zeros(2*(num_el  ements+1),  1); 
whi  I  e  add_mass  ==  1  y '  ; 

mass_node  =  input)'  Node  number  for  lumped  mass  ?  '); 

mass_coord  =  input)1  Translation  or  Rotation  for  lumped  mass  ?  (t/r)  1 

if  mass  coord  ==  1  t 1  ;  %  Translational  lumped  mass 
mass  DOF  =  2  *  mass  node  -  1; 
elseif  mass_coord  ==  'rr;  %  rotational  lumped  mass 
mass  DOF  =  2  *  mass_node; 
end 

mass  d  i  a  g  (  mas  s_  DOF )  =  input)'  Enter  value  of  mass/inertia  (in  "lbf-secA2/ 
%'puts  lumped  mass  on  correct  DOF 
add  mass  =  input)'  Add  another  lumped  mass  ?  (y/n)  ','s'); 

°7o  can  continue  adding  mass  until  1  n '  is  entered 

end;  %  End  while  loop 

end;  %  End  exist)'  ma s s _ d i  a g '  ) 

mx  =  mx  +  d i  a g (  ma s s _ d i  a g ) ;  %  A d d  I  u mp e d  masses  to  [  mx ]  ; 

********************  END  ADDLUMPMASS _j  r  m.  M  *********************** 
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c. 


ASSEMLESENS  JRM 


%  ********************  A  s  s  e  mb  I  e  S  e  n  s  jrm  m  *********************** 

% 

%  This  program  asse mb les  the  total  sensitivity  matrix,  T_ s e n s  tot  and 
%  total  lam  vector,  vect  lam  tot  and  asse mb  les  the  relative  Frequency 
%  error  between  the  natural  frequencies  of  Beam  A  and  Beam  X 
%  Written  by  Constance  R  S  Fernandez,  Spring  2004 
%  Updated  by  John  R  Mentzer,  Spring  2007 


%  I  NPUTS 


%  v  e  c  t  _  I  a  mx  o  s  e  t 
%  v  e  c  t  _  I  a  m  o  s  e  t 
%  v  e  c  t  _  I  a  m" 

%  T  sens  oset 
%  T'sens" 

%  ' 

%  OUTPUTS 

% . 

%  vect  OSET 
%  v  e  c  t  _  I  a  m  tot 
%  T  sens  tot 
%  freq  OSET 
%  f  r  eq"OSETx 
%  rel _f  r  e  q  ERROR 

vect  OSET  =  v  e  c  t  _  I  a  mx  oset  -  v  e  c  t  _  I  am  oset; 

%  I  a  mx  from  actual  beam  with  error  oset,  lam  from  FE  model  oset 
%  Creating  a  vector  of  lam  differences  calculated  ( Lx-  La) 
if  vect  OSET  ==  0; 

%wFien  vector  is  empty  at  first,  the  total  vector  is  equal  to  the 
%  lam  vector  of  Beam  A,  i  .  e .  ,  the  first  19  values  of  vect  I  a  m_  t  o  t  are 
%  the  natural  freq  squared  (radA2/secA2)  of  BeamA 

vect  lam  tot  =  v  e  c  t  _  I  am; 

else 

vect_l  am_tot  =  c  a  t  ( 1 ,  v  e  c  t  _  I  am,  vect  OSET); 
end 


if  T_sens_oset  ==  0; 

T  sens  tot  =  T  sens; 

else 

T  sens  tot  =  c a t ( 1 ,  T  sens,  T  sens  oset); 
end’’ 

freq  OSET  =  s  q  r  t  ( a  bs  ( vec  t  OSET))/2/pi; 

%  Natural  frequency  vector  of  BeamA  in  Hz 
freq  OSETx  =  sqrt(abs(vect_lamx_oset))/2/pi  ; 

%  Natural  frequency  vector  of  Beam  X  in  Hz 
rel  f  r  eqERROR  =  freq  OSET. /freq  OSETx*100; 

%  Relative  error  between  Beam  A  OSET  natural  freq  and  Beam  X  OSET  natural  Freq. 

o/0  ********************  END  Asse  mb  I  e  S  e  n  s  jrm.  m  *********************** 
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BEAMPROPERTIES  JRM 


%  ********************  BeamPropert  i  es  jrm.  m  *********************** 

%  This  is  the  “props  tile"  to  load  nominal  beam  data. 

%  This  program  is  called  by  BeamA  Prompters  to  provide  beam  properties 
%in  order  to  build  BeamA. 

%  This  program  was  written  by  Constance  Fernandez,  Spring  2004 
%  This  progam  modified  by  John  R  Mentzer,  Spring  2007 
%  Out  put  s 

% . 

%  depth 
%  width 
%  E 
%  r  h  o 

%  t  o  t  a  I  length 
%  n  u  m  e  F  e  me  n  t  s 
%  n  o  mf  n  a  I  _  E I 
%  n  o  mi  n  a  I  _  a  r  e  a 
%  nominal  "density 

%  Following  are  actual  measurements  from  experimental  set-up  cantilever 
%  beam. 

depth  =  .504;%  in  z-dir  (inches) 
width  =  1.506;  %  i  n  y-dir  (inches) 

E  =  1 0  e  6 ; 

%E  =  1.  6  5  e  6 ;  %  I  bf  /  s  ec  A2-  i  n  (  1 0  e  6  -  ksi  ) 

%(  1  bf  /  i  n  *2  =  6  8  9  4.  7  6  P  a )  -  >  E(lbf/inA2)  =  ()Pa/  6  8  9  4.7  6 
rho  =0.  1  1  0  4  6  0  9  3  4;  %0.  0  9  8;  %  I  bf/i  nA3 

%T6  temper  alloys  require  a  3  5  -  ksi  tensile  strength,  3  0  -  ksi  yield 
%  strength  and  a  1 0  e  6  -  ksi  elastic  modulus.  Alloy  6  0  6  1  -  T  6  has  1.0 
%  pet  magnesium,  0,6  pet  silicon,  0.3  pet  copper  and  0,2  pet  chromium. 

%  It  has  a  4  5  -  ksi  tensile  strength  and  3  5  -  ksi  yield  strength, 1  The 
%  machinability  of  aluminum  alloys  are  high  (300)  compared  to  titanium 
%  (40).  Aluminum  alloys  can  easily  be  bent  and  provide  easy  loading  and 
%  unloading  of  parts.  Also,  aluminum  is  a  highly  conductive  metal 
%  c  o  mp  a  r  e  d  to  titanium. 


tall  me  a  s  u  r  e  me  n  t 
total  length 
num  elements 
n  o  mf  n  a  I  _  E I 
no  mi  nal  "area 
nomi  nal  "densi  ty 


of  distance 
=  42; 

=  20; 

I  wi  dt  h 


are  in  i  nches 


depth 
rho;  " 


*  depths  / 

*  w  i  d  t  h ;  %  in' 
I  bf/i  n  A3 


12) 

2 


E; 


*********** 


END  BeamProperti  es_j  rm.  m 


E. 


BEAMX  PROMT  JRM 


%  ********************  Rea mX  Prompt  jrm.  m  *********************** 

%  Wr  i  1 1  e  n  By  Prof  Gor  di  s 

%  Inputs  needed 

% . 

%  e  I  e  me  n  t  _  E I 
%  e  I  e  me  n  t  ^  ma  s  s 

%  Out  put  s 

% . 

%  i  Ibis 

%  change  mass,  change  El 
%  new  I  b  F  s 

%  updated  element_mass  in  column  2 
%  updated  el  e  me  n  t  _  E I  in  column  2 
%  ma  s  s  Ibis  ■  index  for  s  e  n  s  i  t  i  v  y  ma  t  r  i  x 
%  E I  Ibis  -  index  for  sensitity  matrix 
%  dv'EI 
%  dv  mass 
%  d  v "  t  o  t 


Prompt  User  for  BeamX  Modification  Data 


di  sp( '  1  ) ;  d i  s p ( '  1  ) ; 

d  |  s  p  ( 1  Modify  nominal  physical  properties  for  second  beam1) 

%  Adjust  mass  values  for  second  beam: 
i .1  bl  s  =0; 

dv  mass  =[  ] ; 
ma  s  s  Ibis  =  [  ] ; 
c  h  a  n  g  e  _  ma  s  s  =  '  n 1  ; 

change_mass  =  input)'  Modify  single/range  element  mass  values  (  y  /  n )  ?  1  ,  1  s 1  ) ; 

%  user  i  nput 

whi  I  e  change_mass  -=  '  n ' ; 

d  i  s  p  ( '  Enter  element  label(s)  for  mass  modification1) 
dispi1  Use  MATLAB  vector  f  o  r  ma  t  >  1  3  5:7  9  ') 

newjbls  =  input)'  >>  '  ,  '  s '  ) ; 

newjbls  =  eval  (['  [' ,  new_l  bl  s, '  ]'  ]);  %  Converts  string  to  vector  of  labels 
i  _  I  b  I  s  =  i  _  I  b  I  s  +  1; 

%  CBS  addition 

mass_l  bl  s(i  _l  bl  s,  1:  I  ength(new_l  bl  s))  =  newjbls;  %  index  for  sensitivity  matrix 

d  i  s  p  ( '  Enter  mass  change  for  element  range') 

mass_change  =  input)'  Enter  percentage  mass  change  ( +/ ■  %)  '); 

dv  mass)  i  Ibls.l)  =  mass  change/100;  %  vector  of  mass  changes  to  second  beam 
( BeamXj 

element  mass(new  I  bis,  2)  =  element  mass(new  I  b  I  s ,  2 )  +.  .  . 

(  mass_change7 100)  *  el  ement_mass(new_l  bfs,  2); 

di  sp( '  '  ) 

change  mass  =  input)'  Modify  another  element  mass  value  (  y  /  n )  7  '  ,  '  s '  ) ; 
di  s p ( 1  '  ) 

end;  %  end  whi I  e 

%  Adjust  El  values  for  second  beam; 
i  J  bl  s  =0; 

change_EI  =  '  n '  ; 
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d  i  s  p  ( '  1  ) 

change_EI  =  input)'  Modify  single/range  element  El  values  ( y  /  n )  ?  ','s'); 
while  change  El  ~=  1  n'  ; 


d  i  s  p  ( 1  Enter  element  label(s)  for  El  modification') 
dispi1  Use  MATLAB  vector  f  o  r  ma  t  >  1  3  5:7  9  ') 

new_lbls  =  input)'  >>  '  ,  '  s'  ) ; 

newjbls  =  eval(['[',new_  ibis,  ']']);  %  Converts  string  to  vector  of  labels 
i  _  I  b  I  s  =  i  _  I  b  I  s  +  1; 

El_l  bl  s(i_l  bl  s,  1:  I  ength(new_l  bl  s))  =  newjbls;  %  index  for  sensitivity  ma  t  r 

d  i  s  p  ( 1  Enter  El  change  for  element  range') 

Exchange  =  input)'  Enter  percentage  El  change  (+/■  %)  '); 

dv_EI  (i_l  bl  s,  1)  =  El  _change/  100;  %  vector  of  El  changes  on  second  Beam 

element  El(newlbls,2)  =  element  El(newlbls,2)+.  .. 

( El  “change/ 100)  *  element_EI[new_lbfs,2); 

d  i  s  p  ( 1  '  ) 

change  El  =  input)'  Modify  another  element  El  value  ( y  /  n )  ?  ','s'); 
di  sp( 1  "'  ) 
end;  %  end  whi I  e 

dv  tot  =  [dv  mass;  d  v  El]; 

%  vector  of  total  changes  to  second  beam  (BeamX)  but  not  location, 

%  End  B  e  a  mX  P  r  o  mp  t .  m 
clear  Exchange  mass_change 


***********  END  Be  a  mX_  P  r  o  mpt  _  j  r  m.  m 
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BOUND ARYCONDITION S  JRM 


o/0  ********************  BoundaryCondi  ti  ons  jrm.  m  *********************** 
%  Wr  i  1 1  en  by  Prof  Gor  di  s 

%  This  script  prompts  the  user  boundary  condition  information 
%  The  script  creates  a  vector  of  DOF  (with  respect  to  the  unrestrained 
%  structure)  and  then  extracts  the  rows  and  columns  of  the  complementary 
%  DOF. 

%  Script  defines  vector  "free_dof  set"  containing 
%  listofunrestraineddof. 

%  The  boundary  conditions  are  applied  in  this  script. 

%  Inputs  needed: 

% . 

%  ndof 

%  k  a ,  ma ,  k  x ,  mx 

%  Out  put  s : 

% . 

%  f  r  e  e  _  d  o  f  set 
%  updated  l<a,  ma ,  k x ,  mx 
%  i  c  n  t  d  o  f 
%  add  3 of 
%  be  node 
%  b  c "  c  o  o  r  d 
%  be' DOF 
%  b  c  boolean 
%  alf  d  o  f  s 
%  restraint  s  wi  t  c  h 

0^  *********7***************************************** 

%  S  t  a  r  t  code: 


if  exist('free_dof_set')==0;  %  Build  f  ree_dof_set  vector 


di  sp( 
di  sp( 
di  sp( 
di  sp( 
di  sp( 
di  sp( 


Sel  ect  a  boundary  condi t i  on  set :  1  ) 

(1)  Cl  amped-free1  ) 

(2)  Cl  amped-CI  amped1  ) 

(3)  Pi  nned-Pi  nned1  ) 

( 4 )  User- Defi  ned1  ) 

( 5)  Free-  Free1  ) 


BC_Choi  ce  =  input)'  >>  Enter  choice:  '); 

if  BC  Choice  ==  1;  %  Cl  amped-free 

free  dof_set  =  [ 3 :  n d o f ] ; 
restrai  nt_swi  tch  =  1  y 1  ; 

el  seif  BC_  Choice  ==  2;  %  Cl  amped-CI  amped 

free  d  o  f  _  s  e  t  =  [ 3:  ndof -  2  ] ; 
restrai  nt_swi  tch  =  1  y 1  ; 

elseif  BC_ Choice  ==  3;  %  Pinned-Pinned 


free  d  o  f  _  s  e  t  =  [ 2 :  ndof-  2  ndof]; 
restrai  nt_swi  tch  =  1  y 1  ; 

elseif  B C_ C h o i  ce  ==  4;  %  User-Defined 


i  c  nt  dof  =0; 

add  dof  =  '  y ' ; 

whi  Te  a d d _ d o f  ==  '  y 1  ; 

b  c  _  n  o  d  e  =  input)'  Node  number  for  restraint  ?  "0"  to  end:  '); 

i  f  b c _  n o d e  ==  0; 
break 

end; 

bc_coord  =  input)1  Translation  or  Rotation  ?  (t/r)  '  ,  '  s '  ) ; 

ient  dof  =  ient  dof  +  1; 
if  bc_coord  =  =  r  t '  ; 

be  DOF) i  c n t  dof)  =  2  *  b c _ n o d e  -  1; 
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elseif  bccoord  ==  '  r 1  ; 

be  DOF  (ient  dof )  =  2  *  b  c  _  n  o  d  e ; 
end;  _  %  End  i  f :  e  I  s  e  block 

end;  %  End  while  a d d _ d o f 

be  boolean  =  ones(ndof,l);  %  [  1  1  1  ...  ient  dof] 

bc'bool  e  a  n  (  be  DOF)  =  zer  os(  I  engt  h(  be  DOF) ,  1 ) ;  %  Put  zeros  in  restrained  dof 
alf_dofs  =  [lTndof];  %  L i  s t  of  all  dof 

free  dof  set  =  all  dofs(l  ogi  cal  (be  boolean));  %  Extract  free  dof 
restrai  nt_swi  tch  =~ '  y 1  ; 

elseif  BC_ Choice  ==5;  %  Free-free  beam 

free  d  o  f  _  s  e  t  =  [1:  ndof  ] ; 
restrai  nt_swi  tch  =  '  n1  ; 

end;  %  End  if-elseif  choice  block 


end;  %  End  exist  block 

ka  =  k  a ( f  r  e  e  dof_set, free  dof  set); 
ma  =  ma(free"dof_set,free~dof~set); 
kx  =  kxjfree'dof  s et , f r ee”dof " s et ) ; 
mx  =  mx(free"dof_set,free"dof"set); 

%  ****************  END  BoundaryConditions  j  r  m.  m  *********************** 
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G. 


BUILD2BEAMS  JRM 


o/0  ********************  Bui  I  d  2  B  e  a  ms  j  r  m.  m  *********************** 

clear 
c  I  c 

%  Revision  history: 

% . 

% 

%  Ver.  1.0:  9/  2  2  /  9  4  Basic  two  beam  assembly 

%  2.0:  Added  multi-element  changes 

%  2.1  3/28/95  Added  read/write  to  file,  rebuild  capability 

%  2.2  3/  2  9/  9  5  Added  lumped  mass  additions 

%  3/10/04  Added  Sensitivity  ma trices,  error  prediction,  plots 

%  2.3  5/10/07  Added  Error  predictor 


% 

%  Program  Description: 

% . 

% 

%  This  program  assembles  the  mass  and  stiffness  ma  trices  for  2  free-free 
%  beams,  referred  to  as  "BeamA"  (analysis)  and  "BeamX"  (experimental).  The 
%  program  can  be  run  in  several  modes: 

% 

%  "Build"  mo  d  e : 

%  The  user  provides  baseline  data  for  BeamA,  assumed  to  be  a 
%  homogeneous,  uniform  beam.  Data  provided: 

% 

%  (1)  Beamlength 

%  ( 2 )  N u mb e r  of  e I  e me n t  s 

%  (  3 )  N  o  mi  n  a  I  El 

%  (4)  Nominal  cross-sectional  area 

%  ( 5 )  No  mi  n  a  I  we  i  g  h  t  density 

% 

%  The  program  then  prompts  the  user  for  instructions  on  how  to  modify 
%  "BeamA"  data  to  arrive  at  "BeamX"  data,  The  user  can  modify  element 
%  masses,  and/or  element  El  values.  The  modification  can  be  applied  to 
%  either  a  single  element,  or  range  of  elements,  e.g. 

% 

%  Modify  single/range  element  mass  values  (  y  /  n )  ?  y 

% 

%  If  "y"  is  entered,  the  user  enters  the  number  of  the  element  for  mass 
%  a  d  j  u  s  t  me  n  t : 

% 

%  Enter  element  label(s)  for  mass  modification:  1 

%  Use  MAT  LAB  vector  f  o  r  ma  t  >  1  3  5:7  9 

% 

%  Enter  percentage  mass  change  ( +/ -  %) 

% 

%  The  user  is  pro mp ted  to  modify  another  element  or  range  of  elements: 

% 

%  Modify  another  element  mass  value  (y/n)7  y 

% 

%  This  process  continues  until  the  user  enters  an  "n"  for  no  change. 

%  This  entire  process  can  then  be  repeated  for  El  adjust  me  nt. 

% 

%  The  program  saves  the  beam  definition  data  in  a  binary  (.mat)  file 
%  "beamdata"  at  the  end  of  execution, 

% 

%  The  program  can  also  be  run  in  "Read"  mode  by  entering  an  " r '  at 
%  the  initial  p  r  o  mp  t , 

% 

% 

%  Script  Execution  Path: 

% 

% 

% 

%  Bui  I  d2Beams_j  rm.  m  --  User  executes  this  program. 

%  BeamA  Prompt  jrm,  m  --  Prompts  User  for  BeamA  nominal  beam  data 

%  Bea  mX~  P r  o  mpt  ~ j  r  m,  m  --  Prompts  User  for  BeamX  modification  beamdata 

%  AssemBl  e2Beams  jrm,  m  --  Called  by  Bui  I  d2Beams,  builds  [ka]  [ma]  [kx] 

%  [  mx ] ,  p  I  o  t  s  f  r  e  q  s , 

%  AddLumpmass  jrm.  m  --  Prompts  User  for  BeamX  lumped  mass  addition 

%  BoundaryConditions  jrm.  m  --  Prompts  user  for  B.C.'s  and  applies  them. 

%  PI  ot  Bea  mModes  jrm,  m  --  Calculate  beam  modes  and  plot  frequencies 
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% 

%  BeamSensi  ti  vi  ty  j  r  m.  m 

%  BeamSensi  ti  vi  tyDSET  j  r  m.  m  - - 

%  recorded  H  jrm.  m 

applied 

%  A  s  s  e  mb  I  e  S  e  n  s  j  r  m.  m 

errors. 

%  ABCr  unTHRU.  m 

%  Saves  data  to  "  b  e  a  md  a  t  a .  ma  t  ‘ 

% 

%  Start  code: 

%  . 

0^  ************************************************************************* 


Calculate  sensitivity  matrix  T-sens 
Calculate  sensitivity  matrix  using  ABC 
Calculates  the  nat.  freq  of  BeamX  with  ABC 

Assembles  the  sens  matrices  and  calculates 

Calculates  the  DV  and  cond  number  of  matrix  used 


d  i  s  p  (  1  Building  2  beams  from  scratch, 


Be  a  mA_  Pr  ompt  jrm; 

Be  a  mX~  Pr  o  mpt "  j  r  m; 

A  s  s  e  mb  I  e  2  B  e  a  ms  _  j  r  m; 
AddLumpmass  jrm; 

% 


%  Prompt  for  BeamA  Data:  run  prompt  script 

%  Prompt  for  BeamX  Modification  Data: 

%  Run  script  to  a  s  s  e  mb  I  e  mass  and  stiffness  ma  t  r  i  c  e  s 

%  B e a mX  I  u mp e d  mass  vector  construction  and 

application 


kx  beam  =  kx;  %  saves  the  BeamX  matrices  without  BC  to  be  used  later 
mx  "  b  e  a  m  =  mx ; 


Boundar  yCondi  t  i  o n s  _ j  r  m; 


Prompt  for, 


and  apply  boundary  conditions 


kx 

mx 

beamBC  =  kx; 
beamBC  =  mx; 

%  saves 

t  he 

Beam 

X 

ma  t  r  i  c  e  s 

wi  t  h 

BC 

t  0 

be 

used 

later 

k  a" 
ma" 

beamBC  =  ka; 
beamBC  =  ma; 

%  saves 

t  he 

Beam 

A 

ma  t  r  i  c  e  s 

wi  t  h 

BC 

t  0 

be 

used 

later 

P I  ot  Be  a mModes  c  r  s 


%  Calculate  beam  modes  and  plot  frequencies 


BeamSensitivity  jrm;  % 
BeamSensi  ti  vi  tyOSETJ  rm;  % 


Calculate  sensitivity  matrix  T-sens 
Calculate  sensitivity  matrix  using  ABC 


recorded  H  jrm; 
Assembl  eSens  jrm; 
ABCr  unTHRU  jrm; 
pi  ot  t  i  n g  BARS  j  r  m; 
error  predictor; 


%  Calulates  the  nat.  freq  of  BeamX  with  ABC  applied 
%  Assembles  the  sens  matrices  and  calculates  errors. 
%  Calculates  the  DV  and  cond  number  of  matrix  used 
%  Bar  plots  of  predicted  DV  vs.  true  error 
%  Generates  error  prediction  for  desired  scenario. 


%  Save  Defining  Parameters  for  Beams  and  plots 

di  s  p( '  ...saving  beam  data  to  file1) 
save  beamdata.  mat 


d i  s p ( 1  Bui  I  d2Beams  end,1) 
% 


o/0  ********************  END  Bui  I  d2Beams  jrm.  m  *********************** 
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H.  DISPLACEMENTPLOT  OSET  JRM 


o/0  ********************  d  j  s  p  |  acementPI  ot  OSET  j  r  m.  m  *********************** 

%  Wr  i  1 1  e  n  by  Constance  Fernandez  Spring  2004 
%  Updated  by  John  Mentzer,  Spring  2007 

%  This  program  plots  the  mode  shapes  (phi,  lam)  phi  vs  nodal  position 
%  of  beam  when  ABC  are  applied, 

%  Inputs 

% . 

%  p  I  o  t  k  x 
%  k  a  0  base 
%  n  u  m  mo  d  e  s  0 
%  phi  moT 
%  phi  APLOT 
%  pinned 
%  n  u  m_  e  I  e  me  n  t  s 

%  Out  put  s 

% . 

%d  i  s  p  1 ,  d  i  s  p  I  a 

%i  i .  g 

%  ypos 
% . 

%  d  i  s  p  1  =  zeros(cei  I  (.  5*si  ze(pl  otkx,  1)),  num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector, 
di  s  pi  =  zeros(cei  I  ((20/38)*size(pl  otkx,  1)),  num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector. 

%  d  i  s  p  1  a  =zeros(ceil(.5*size(kaO  base,l)),num  modesO); 

%i  n  i  1 1  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector, 
displa  =  zeros(  cei  I  ( ( 20/38)*si  ze(  kaO  base, 1)), num  modesO); 

%initilize  disp  vector  and  provides  the  first  zero  of  the  vector, 
for  jj  =  1:  cei  l((20/38)*si  ze(pl  otkx,  1)); 

di  s  p  1  { j  j  +1,  :  j  =  phi  XPLOTj  2*j  j  ■  1,  1:  num_  modesO) ; 

%  every  other  phi  to  give  displace  me  nt  at  sequential  nodes 
di  spla(jj  +1,  :  )  =  phi  A  P  L  OT  (  2*j  j  -  1,  1;  num  modesO); 
end 

%  This  loop  normalizes  the  modes  shapes  to  the  tip  mo  dal  displace  me  nt. 
if  pinned  ==  41  %  tip  pinned  is  a  special  case,  no  new  calculations  are  needed 
di  s  p 1 (  :  ,  :  )  =  di  s  p  1  (  :  ,  :  ) ; 
d  i  s  p  1  a  ( :  ,  ;  )  =  displa(;,;j; 

elseif  I  engt  h(  pi  nned)  >1  &  pi  nned(  1,  2)  ==41  %t  i  p  pinned  with  2  ABC's 
d  i  s  p  1  (  :  ,  :  )  =  d  i  s  p  1  (  :  ,  :  ) ; 
d  i  s  p  1  a  (  :  ,  :  )  =  displa(:,:); 
else 

for  g  =  l:num  modesO 

d  i  s  p  1  ( :  ,  g )  =  di  spl(  ;  ,  g)/di  spl(num  el  eme  nt  s  +1 ,  g ) ; 
d  i  s  p  1  a  ( :  ,  g )  =  di  spla(  :  ,  g)/di  spla(num_el  ements+1,  g); 
end 
end 
%  end 

%y  p  o  s  =  [  1 :  1 :  num  el  ements+1] ;  %  Location  of  nodes  used  in  plotting 
ypos  =  [  1 :  1 ;  n  u  m  el  ement  s+1] ;  %  Location  of  nodes  used  in  plotting 
%  i  f  mass  Ibis  ~=  [  ] ; 

% 

%  for  k  k  =  1 ;  s  i  z  e  (  ma  s  s  I  b  I  s ,  1 ) ; 

%  ff  =0; 

%  for  JJ  =  1:  I  engt  h(  f  i  n  d  (  mass  I  bl  s(  kk,  :  )  >0) ) ; 

%  ff  =  ff+1; 

%  p  o  s  m(  k  k ,  2  *  J  J  - 1 )  =  ma  s  s  I  b  I  s  ( k  k ,  f  f ) ; 

%  posmjkk,  2  *  J  J  )  =  mass  I  51  s(kk,  ff)  +1; 

%  end 

%  end 

% 

%  i f  kk  ==  1 
%  p  o  s  M  =  p  o  s  m; 

% 

%  else 

% 

%  f  o  r  u  u  =  1 :  k  k  -  1 ; 

%  posM  =  cat(2,  p  o  s  m(  u  u ,  :  )  ,  posm(  uu+1,  :  )  )  ; 

%  end 

% 

%  end 

%  posM  =  s  or  t  ( pos  M(  f  i  n  d  (  pos  M>0) ) ) ; 
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%  m  =  .  5 * o n e s ( si  ze(  posM) )  +i  cnt  oset; 

%  end 
% 

%  i  f  El  Ibis  ~=  [  ] ; 

%  for  kk  =  1 :  si  z  e  (  El  I  b  I  s ,  1 ) ; 

%  ff  =0; 

%  f or  J  j  =  1:  I  engt h(f  i  nd(  El  I  bl  s(  kk,  :  )  >0) ) ; 

%  ff  =  f  f  +1; 

%  p  o  s  e  ( k  k ,  2  *  j  J  - 1 )  =  E I  I  b  I  s  ( k  k ,  f  f ) ; 

%  p o s e ( k k ,  2* J  J  )  =  El  I  Bl  s(  kk,  f f )  +1; 

%  end 

%  end 

% 

%  i  f  kk  ==  1 
%  posE=  pose; 

% 

%  else 

% 

%  for  uu  =  l:kk-l; 

%  posE  =  cat(2,  p  o  s  e  (  u  u ,  :  )  ,  pose(  uu+1,  :  )  )  ; 

%  end 

%  end 

%  posE  =  sort(posE(find(posE>0))); 

%  e  =  - .  5*ones(  si  ze(  posE) )  +i  cnt_oset ; 

%  e  n  d 


o/0  ********************  END  displace  me  n  t  P I  o  t  OSET  jrm.  m  *********************** 
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ERROR  PREDICTOR  JRM 


%**************g|-|-Qp  predictor  jrm.  m* *************** 

%Mode  error  predictor 

%Wr  i  1 1  e  n  by  John  Mentzer,  Spring  2007 

%  Dv=zeros(  40,  20)  ; 

%  f  o  r  n  =  1 :  4  0 ; 

%  Dv(n,:)=(T  sens  tot(n,  :  )'  )*((l  amx(n,  :  )/(2*pi  ))-(l  ama(  n,  :  )/(2*pi  ))); 

%  end 

f  o  r  ma  t  long 

%///////////  ORGI  Nl  AL  METHOD///////////////////// 
for  q  = 1 :  2  0 

dv_cal  cul  ated  b  a  s  e  2  (  :  ,  q )  =  T_  s  e  n  s  t  ot  ( 1:  q,  :  )  \  vect  I  am_tot(  1:  q); 
end 

for  t  =40:  59 

dv_cal  cul  ated_base3(:  ,  (t-39))  =  T  sens  tot(40:t,:)\vect  lam  tot(40:t) 
end 


%f  o  r  BASE 
for  p  = 1 : 2  0 

I  i  mi  t_cond(  1,  p)  =cond(T  sens  tot(l:p,:)j; 
end 

%f  o  r  ABC  system 
for  AB Cp  =4 0 :  59 

I  i  mi  t_condABC(  1,  (ABCp  -  39)  )  =cond(T_sens  tot(40:ABCp,  :  )); 
end 


%  for  g  =1 :  20 

%  dv  calculated  base_trans(g,  :  )  =  vect  lam  tot(g,  :  )*(T_sens_tot(g,  :  )) 
%  end 


%  f  i  g  u  r  e  ( 8 ) 

%  w=l:  2  0 

%  b  a  r  (  w,  v  e  c  t  I  am_tot(w,  1),  .  25,  '  r1  ) 
%  %  % 


d  v  c  a  I  c  u 
dv'cal  cu 
dv'cal  c  u 
dv'cal  c  u 
dv'cal  cu 


at  ed_basell( 
at  ed_base21( 
at  ed_base31( 
at  ed_base41( 
at  ed_base51( 


,  1) 
,  2) 
,  3) 
,  4) 
,  5) 


=  T_ s  e  n s  t  ot (  1, 
=  T_sens"t  ot  j 2, 
=  T_sens~t  ot ( 3, 
=  T  sens“tot(4 


)  \  vect  I  am  t  ot  ( 1) 
)\vect"lam"tot(2) 
)  \  v  e  c  t  _  I  a  m_  t  o  t  ( 3 ) 
)  \  v  e  c  t '  I  a  m_  t  o  t  (  4 ) 


=  T_sens“t  ot  (  1:  5,  :  )  \  v  e  c  t  _  I  a  m_  t  o  t  ( 1 :  5); 


d  v  _  c  a  I  cul  at  ed_base2  /  ( norm(  dv_cal  cul  at  ed_base2(  :  ,  :  ),  i  n  f ) ) ; 


o^*******************Efld  of  error  predictor  jrm.  m* ******************* 
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J.  FBEAMKM  JRM 


%  %  ********************  f  be  a  mk  m  j  rm  rn  *********************** 

% 

%  %  function  [  k  b  e  a  m,  mb  e  a  m]  =f  b  e  a  mk  m(  I  ,  e  i  ,  m) 

%  %  Provided  by  Prof  Gordis 
% 

%  function  [  k  b  e  a  m,  mb  e  a  m]  =f  b  e  a  mk  m(  I  ,  e  i  ,  m) 

%  % 

%  % 

%  %  This  function  returns  the  stiffness  and  mass  ma trices  for 
%  %  a  si  mp I  e  2 -  n o d e  b e a m  e I  e me n t . 

%  % 

%  %  Note:  m=rho  *  area  *  length  =  total  element  mass 
%  % 

%  %  Reference:  R.  D.  Cook,  Concepts  and  Applications  of  F.E.  Analysis 
% 

%  %  Out  put  s 
%  % . 

%  %  k  b  e  a  m,  mb  e  a  m,  i  ,  j 
% 


%  k  be  a  m=z  e  r  os  (  4,  4 ) ; 

%  mbe a m=z  e r  os  (  4,  4 ) ; 

%  % 

%  kbeam(  1,  1)  =12.  0; 

%  k b e a rm(  1,  2)  =6.  0*1  ; 

%  k b e a m(  1,  3)  =-  12.  0; 

%  kbeamj  1,  4)  =6.  0*1  ; 

%  k b e a m(  2,  2)  =4.  0*1  A2; 

%  k b e a m(  2,  3)  =■  6.  0*1  ; 

%  k b e a m(  2,  4)  =2.  0*1  A2; 

%  k b e a m(  3,  3)  =12.  0; 

%  k b e a m(  3,  4)  =■  6.  0*1  ; 

%  kbeamj  4,  4)  =4,  0*1  A2; 

%  % 

%  mb e a m(  1,  1)  =156.  0; 

%  mb e a m(  1,  2)  =22.  0*1  ; 

%  mbe amj  1,  3)  =54.  0; 

%  mbeam(  1,  4)  =-  13.  0*1  ; 

%  mbe  a  mi  2 ,  2  j  =4 ,  0*  I  A2 ; 

%  mb  e  a  m(  2 ,  3)  =13.  0*1  ; 

%  mbeam(  2,  4)  =-  3.  0*1  A2; 

%  mbeamj  3,  3)  =156.  0; 

%  mb  earn)  3,  4)  =-  22.  0*1  ; 

%  mbeamj  4,  4)  =4.  0*1  A2; 

%  % 

%  f  o  r  i  =  1 :  4 ; 

%  for  j  =i  :  4 ; 

%  k  be  a  m(  j  ,  i  )  =k  bea  m(  i  ,  j  ) ; 

%  mb  earn)  j  ,  i  )  =mbeamj  i  ,  j  ) ; 

%  end 

%  end 
%  % 

%  k b e a m=(  ei  / 1  A3 )  *  kbeam; 

%  mb  e  a  m=  i  m/  4  2  0 ,  0  )  *  mb  e  a  m; 

%  % 

%  %  e  n  d  function  b  e  a  mk  m 
% 

%  %********************  END  fbeamkm  jrm.  m  *********************** 
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K.  FMODES 


%  %  ********************  f  Mod es  m  *********************** 

% 

%  function  [lam,  phi  ]=fmodes(k,  m,  n  u  m  to  print); 

%  %  Provided  by  Prof  G o  r  d  i  s 

%  %  This  program  prints  to  the  screen  natural  modes  of  system  (phi). 

%  % 

%  %  Usage;  [lam,  phi]=fmodes(k,m,  num  to  print) 

%  %  '  ' 

%%  This  function  can  be  used  with  1  to  3  arguments,  as  follows; 

%  % 

%%  [  I  am,  phi  )  =f  modes!  a)  :  Produces  modes  of  [a]  with  no  print  of  freqs  in 

Hz . 

%  %  [I  am,  phi  ]  =f  mo  d  e  s  ( a ,  i  )  :  Produces  mo  d  e  s  of  [  a  ]  with  print  of  "  i  "  freqs  in 

Hz . 

%  %  [I  am,  phi  ]  =f  mo  d  e  s  ( k ,  m)  ;  Produces  modes  of  [m\k]  with  no  print  of  freqs  in 

Hz. 

%  % 

%  % 

%  % 

%%  This  function  returns  a  vector  containing  eigenvalues  (rad/ sec)  A2 
%  %  and  a  matrix  containing  the  mass  normalized  mode  shapes, 

%  %  The  mode  information  is  sorted  by  frequency  in  ascending  order. 

%  %  If  num_to  print  >  0;  tabular  listing  of  numto  print  freqs  in  Hz  is 

printed. 

%  %  If  numto  print  <=0,  no  print, 

%  '  ' 

%  %  Inputs 

%  % . 

%  %  v,  i  ndex,  m,  k 
% 

%  %  P  r  o  g  r  a  ms 

%  % . 

%  %  f  N  o  r  ma I  i  z  e 
% 

%  %  Out  put  s 

%  % . 

%  %  phi 

%  %  n  u  m_  t  o  print 
%  %  e  r  r  o  r 
%  %  e 
% 

%  % 


% 

%  i  f  n  a  r  g  i  n  ==  1; 

%  [A]  w /  no  print  request  for  freqs  in  Hz. 

% 

%  v(  1,  :  )  =  1  normal  i  zati  on 
%  ( v,  d]  =ei  g(  k) ; 

%  [temp,  i  nai  ces]  =  s  o  r  t  ( a  b  s  ( d  i  a  g  ( d ) ) ) ; 

%  lam  =  diag(d); 

%  lam  =  lam(indices); 

%  [  phi  ]  =f  Nor  mal  i  z  e(  v( :  ,  i  ndi  ces ) ,  'one'); 

%  num  to  print  =  0; 

%  '  ‘ 

%  elseif  nargin  ==  2  &  size(m,l)  ==  1;  %  [A]  w/  print 

request  for  freqs  in  Hz. 

% 

%  v(  1,  :  )  =  1  n o r  ma I  i  z  a t  i  o n 
%  [  v,  d]  =ei  g(  k) ; 

%  [  temp,  i  ndi  ces]  =  s  o  r  t  ( a  b  s  ( d  i  a  g  ( d ) ) ) ; 

%  lam  =  diag(d); 

%  lam  =  lam(indices); 

%  [phi]=fNormalize(v(:,  indices),  'one'); 

%  n  u  m  t  o  p  r  i  n  t  =  m; 

%  '  " 

%  e  I  s  e  i  f  n  a  r  g  i  n  =  =  2  &  s  i  z  e  (  m,  1 )  >  1 ;  %  [  k  ] ,  [  m] 

w /  no  print  request  for  freqs  in  Hz. 

% 

%  mass  n  o  r  ma  I  i  z  a  t  i  o  n 
%  [  v,  d]  =ei  g(  m\  k) ; 

%  [l  am,  i  ndex]=sort(abs(di  ag(d))) ; 

%  [  phi  ]=fNormal  i  ze(v(:  ,  i  ndex),  1  mass'  ,  m); 

%  num  to  print  =  0; 

%  '  “ 
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[  k] ,  [  m]  w/  print 


%  e I  s e i  f  n a  r  g i  n  ==  3  &  s  i  z  e ( k ,  1 )  >  1  &  s  i  z  e (  m,  1 )  >  1 ;  % 
request  for  freqs  in  Hz, 

% 

%  mass  normal  i  zati  on 
%  [  v,  d ]  =ei  g(  m\  k) ; 

%  [I  am,  i  ndex]=sort(abs(di  ag(d))); 

%  [  phi  ]=fNormal  i  ze(v(:  ,  i  ndex),  1  mass'  ,  m); 

% 

%  else 
% 

%  num  to  print  =  - 1 ; 

%  error(rError  from  f  Modes,  m:  Check  input  arguments.') 

% 

%  end 
% 

%  i  f  n  u  m  t  o  print  >  I  e  n  g  t  h  ( k )  ; 

%  _num  to_print  =  length(k); 

%  end 
% 

%  if  nargin  <  3  &  rem(  I  engt  h(  k) ,  2)  ==0  &  k(  1:  I  engt  h(  k)  /  2,  1:  I  engt  h(  k)  /  2)  == 
zeros(length(k)/2,length(k)/2);  %  Have  [A]  matrix 
%  e  =  1;  %  Ei  genval  ues  are  wn 

%  else 

%  e  =  0 .  5 ;  %  Eigenvalues  are  wnA2 

%  end 


%if  num  to  print  >0; 

% 

%  d  i  s  p  ( '  1  ) ,  d  i  s  p  ( 1  '  ) 

%  d  i  s  p( ' - 1  ) 

%  di  sp( '  Freqs  in  Hz .  :  1  ) 

%  d  i  s  p  ( ( I  a  m(  1 :  n  u  m  t  o _  p  r  i  n  t ) .  *e )  /  2  /  p  i  ) 

%  di  s p( 1 - ) 

% 

%  end 
% 

%  %  ********************  END  f  Modes .  m  *********************** 
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L.  FNORMALIZE 


********************  fNormalize.m  *********************** 


function  [phi]  =  fNormal  i  ze(phi  ,  method,  m); 

% 

%  Usage:  [phi]  =  fNormal  i  ze(phi  ,  method,  m); 

% 

%  phi:  matrix  whose  columns  are  to  be  (independently)  normalized 
%  method:  String  variable.  The  following  choices  are  available: 


%  % 


%  % 


mass 

inf' 

one' 

I  e  n  g  t  h 1 


Mass  normal  i  zati  on 
Infinity  normalization 
First  element  =1 
Length  =  1 


m:  matrix  used  for  normalization,  i  .  e .  ,  phi1*  m*  phi  =  eye 


% 

% 

%  % 
% 

% 

% 

%  % 
% 

% 

p  h  i 
% 

% 

% 


phi 


s  wi  t  c h  met  hod 

case  'mass' 


Mass  normal  i  zati  on 


d  I  s  p  ( '  mass  normalization1) 

phi  =  phi  *  diag(sqrt(diag((phi'  * 


case  ' i  nf 1 


m  *  phi  ) .  A(  - 
Infinity  normalization 


dispj'inf  normalization1) 
for  icnt  cols  =  1:  s  i  z  e(  phi  ,  2) ; 
p"h  i  (  :  ,  i  c  n  t  cols)  = 

:  ,  i  c n t  col  s)  /  norm(  phi  ( :  ,  i  cnt_col  s) ,  i  nf) ; 
end 

case  'one'  %  First  element  =  1 

d  i  s  p  ( '  one  normalization1) 
phi  =  phi  *  di  ag(  ( phi  (  1,  :  ) .  A(  -  1) ) 1  ) ; 

case  'length'  %  Length  =  1 

dispj' length  normalization1) 
for  icnt  cols  =  l:size(phi,2); 
pli  i  (  :  ,  i  c  n  t  cols)  = 

:  ,  i  c  n  t  cols),  /  nor  m(phi(:,icnt_cols),  'fro'); 
end 

end 


********************  end  fNormal  i  ze.  m  *********************** 


l)))); 
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M.  FOSET  FROM  ASET 


%  %  ********************  fOset  from  As  6 1  m  *********************** 

%  ■ 

%  function  [oset]  =  fOset  from  Aset(ndof,aset); 

%  % 

%  %  Usage:  [oset]  =  fOset  from  Aset  ( ndof ,  aset ) ; 

%  % 

%  %  This  function  d  e  t  e  r  mi  n  e  s  the  compl  ementary  subset  "oset" 

%  %  from  a  set  [ 1 : 1 :  ndof]  and  the  subset  aset  =  [x  x  x  ...]. 

%  % 

%%  ndof:  Total  number  of  DOF,  Set  is  labeled  "nset". 

%%  aset:  Retained  DOF  (proper  subset  of  [  1 : 1 :  ndof]) 

%  %  oset:  aset  U  o  s  e  t  =  n 
%  % 

%  %  Pr  ovi  ded  by  Prof  Gor  di  s 
%  % 

%  . 

%  n  s  e  t  =  [  1 :  n  d  o  f  ] ; 

% 

%  for  icnt  =  1  :  length(aset); 

%  i  ndi  c  e  s  (  i  cnt )  =  f  i  n  d  (  n  s  e  t  ==  aset(icnt)); 

%  end 
% 

%  bool  =  o  n  e  s ( s i  z  e ( n  s  e  t ) ) ; 

%  bool  (i  ndi  ces)  =  zeros(size(indices)); 

%  oset  =  nset  ( f  i  nd( bool  >0) ) ; 

% 

%  %  ********************  fOset  from  Aset.m  *********************** 
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N.  FSPRINGMASS2 


%  %  ********************  f  Spr i  ngMass  2  m  *********************** 

%  function  [  k,  m]  =fSpri  ngMass2(  spri  ngs,  mass,  BC)  ; 

%  % 

%%  Usage:  function  [k,  m]  =fSpri  ngMass2(spri  ngs,  mass,  BC) 

%  % 

%  %  This  function  script  a  s  s  e  mb  I  e  s  the  stiffness  [  k  ]  and  mass 

%  %  [  m]  ma  t  r  i  c  e  s  for  an  a  s  s  e  mb  I  a  g  e  of  springs. 

%  % 

%  % 

%  %  A  linear  chain  of  springs  and  masses  is  a  s  s  u  me  d . 

%  %  The  number  of  springs  is  defined  by  the  length  of  the  vector 

1  spri  ngs 1 

%  %  and  their  values  by  the  elements  of  'springs.1 

%  % 

%  %  The  number  of  masses  is  defined  by  the  length  of  the  vector  'mass' 

%  %  and  their  values  by  the  elements  of  'mass'. 

%  %  NOTE:  The  number  of  masses  must  equal  to  the  final  number  of 

act i  ve 

%%  DOF  ( i  .  e.  ,  af  t  er  BC' s  a ppl  i  ed) . 

%  % 

%  %  Boundary  conditions  are  specified  by  the  vector  'BC.'  This  vector 
%%  contains  the  DOF  numbers  which  are  to  be  restrained. 

%  % 

%  %  For  example,  to  build  the  following  system: 

%  % 

%  %  .01  .02  .015 

%  %  |  --////--[  m]  --////--[  m]  --////-■[  m] 

%  %  5  6  3.4 

%  % 

%  s  p  r  i  n  g  s  =  [  1  1  ] ; 

%  ma  s  s  =[.01.01]; 

%  BC  =  [  1] 

%  % 


BEGI  N  SCRI  PT 


%  if  length(mass)  ==  ( I  en  gt  h  ( s  pr  i  ngs ) +1 )  -  length(BC); 

% 

%  k  =  z er os  (  I  engt h(  s pr i  ngs )  +1,  I  engt h ( s  pr  i  ngs ) +1) ; 

%  m  =  z  e  r  o  s  ( I  e  n  g  t  h  (  ma  s  s ) ) ; 

% 

%  %  a  s  s  e  mb  I  e  stiffness  ma  t  r  i  x : 

% 

%  r  ows  =  [  0  1] ; 

%  for  i s  p  r i  n  g  =  1  :  I  e  n  g  t  h ( s  p  r i  n  g  s ) ; 

% 

%  rows  =  rows  +  1; 

% 

%  addthis  =  [springs(ispring)  -spri  ngs(i  spri  ng);  -  spri  ngs(i  spri  ng) 

s  pr  i  ngs ( i s  pr i  ng)  ] ; 

%  k  ( rows,  rows)  =  k(  r  ows ,  r  ows )  +  addthis; 

%  end 

% 

%  if  ~i  sempt  y(  BC) ; 

%  keep  =  fOset_from  Aset(  I  ength(  spri  ngs)  +1,  BC)  ; 

%  k=k(keep,keep); 

%  end 

% 

%  %  a  s  s  e  mb  I  e  ma  s  s  ma  t  r  i  x : 

% 

%  m  =  d  i  a  g  (  ma  s  s ) ; 

% 

%  else 

% 

%  d  i  s  p  ( 1  E  r  r  o  r  in  f  S  p  r  i  n  g  ma  s  s  2 .  Check  #  masses,  springs,  and  B  C "  s .  1  ) 

%  return 

% 

%  end 

%  %  ********************  end  f  Spr  i  n g Ma s s 2  m  *********************** 

% 
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O.  NORMRUNTHRU  CRS 


o/0  ********************  normRUNthru  c  r  s .  m  *********************** 

%  This  program  finds  the  NORM  of  the  columns  of  sensitivity  matrix  and  the 
%  NORM  of  the  rows  of  the  inverse  of  the  sensitivity  matrix.  This  was  used 
%  t  o  find  a  correlation  of  the  good  prediction  to  the  ABC  system  used.  It 
%  also  plots  the  information  in  helpful  graphes, 

% 

%  This  program  was  written  for  a  system  of  19  natural  freq.  Set  of  5  modes 
%were  used  in  each  ABC  system,  i  .  e .  ,  ,  modes  1-5,  modes  6  - 10,  or  modes  11-15. 
%  This  accounts  for  the  3  sets  of  modes  per  condition  as  listed  below  in 
%  the  "for"  loop  modeN  =  1:3.  This  program  comp  ares  the  use  of  the  first  5 
%  modes  of  the  base  system  and  one  set  of  five  modes  of  the  ABC  to  that  of 
%  o f  just  10  modes  of  the  ABC.  Notice  that  the  sensitivity  is  a  10x10 
%  square  matrix  indicating  that  only  mass  or  El  changes,  not  both  were 
%  ma  d  e . 

%  This  program  is  not  part  of  Bui  I  d2Beams_crs.  m  program.  It  is  run 
%  s  e  p  a  r  a  t  e  I  y . 

%  Written  by  Constance  R  S  Fernandez,  Spring  2004 

%  Inputs 

% . 

%  cond_basePI  us 
%  i  c  n  t  _  o  s  e  t 
%  T  s  e  n  s  _  t  o  t 
%  E F  Ibis,  mass  Ibis 
%dv~mass,  dv_ET,  dv_cal  _ABCt  en 

%  Out  put  s 

% . 

%  BASE 
%  BASET 

%abcN,  count  N ,  a  cN 
%  modeN,  st  art  modeN,  st  art  modeNT 
%  bb,  t,  tl  NV,  cc,  T,  T I  NV,  vv,  tt 
%  model  a  be  I  NORM 

%  norm  vectT,  norm  vecTABC,  norm_vecTi  nvABC 
%  normC 

%  baseN,  baseABCN,  abc  conN,  abc_conTN 
%  baseABCNt  en,  abc  conRten,  abc_conTNten 


%  ----Start  program - 

BASE  =  int2str(cond  basePI  us(  1) ) ;  %  cond  no.  of  the  base  line  system 
BASET  =  spri  ntf(  1  Basel  1:  5]  Cond  =  %s',  BASE);  %  used  for  plotting 

%  ======|  nitializatio  n  =====% 

abcN  =  0; 
count  N  =0; 

%  ====  =  ==Cal cul  ati  ons  of  NORM  vectors  ====  =  =% 


for  count  =  1 :  i  c  n  t  oset  +1  %  number  of  conditions  (base  +  ABC) 
a  _  c  N  =  1; 

for  modeN  =  1:3  %  3  sets  of  modes  per  boundry  condition 

start  mo  d  e  N  =  abcN  +  a  _  c  N ;  %  the  beginning  mode  number  of  each  set 

%  indicates  the  use  of  the  first  5  modes  of  base  system  +  5  modes 
%  of  ABC  system 

bb  =  [1:5,  start  modeN:  start  modeN +4]; 

%  labeling  of  modes  for  plotting 
model  abel  NORM  =  i  nt2str(a_cN:  a_cN+4); 

a  _  c  N  =  a  _  c  N  +  5 ;  “/.advances  to  the  next  set  of  modes 

%  base  system  plus  5  modes  of  ABC 

t  =  T  sens  t  o  t  (  b  b ,  :  )  ;  %  10x10  matrix 

1 1  NV  =  i  nv[T_sens_t  ot  ( bb,  :  ) ) ;  %  10x10  matrix 

%  f  i  r  s  t  10  modes  of  ABC  solo 

start  mode  NT  =  abcN+1;  %  the  beginning  mode  number  of  each  set 
cc  =  [  s  t  a  r  t  mo de  NT:  s  t  a r  t  mod e NT +9 ] ; 
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T  =  T_  s  e  n  s  t  o  t  ( c  c ,  :  )  ;  %  10x10  matrix 
T I  N  V  =  i  nv[T_sens_t  ot  ( cc,  :  ) ) ;  %  10x10  matrix 

%  for  loop  for  NORM  of  columns  and  rows  of  invfsens  matrix) 
for  vv  =  1:10  %  10  rows,  10  columns 
%  base  +  ABC  system 

nor  m_vecT(  vv,  count  N  +  modeN)  =  nor  m(  t  (  :  ,  vv) ) ;  %  columns 
nor  m_vecTi  n  v  ( vv,  count  N  +  modeN)  =  nor  m(  1 1  N  V  ( vv,  :  ) ) ;  %  rows 
%  f  o  r  10  mo  d  e  s  ABC  solo 

norm_vecTABC(vv,  countN  +  modeN)  =  n  o  r  m(  T  ( :  ,  v  v ) ) ;  %  columns 
nor  m_vecTi  nvABC(  vv,  count  N+modeN)  =  nor  m(  Tl  NV(  vv,  :  ) ) ;  %  rows 

end  %  vv  loop 
end  %  ModeN  loop 

abcN  =  abcN+19;  %  advances  to  the  next  ABC  system 
count  N  =  count  N  +  3:  %  counts  up  each  set  of  ABC 
end  %  count  loop 
nor  mC=4;  %  initialize  for  plots 

%======= ===P L OTT I  NG==========% 

for  1 1  =1 :  10  %figures  (  3  0  -  4  0  )  plots  6  graphes  per  figure 
f i  g  u  r  e ( 1 1  +3  0 ) 
s  u  b  p  I  o  t  ( 3 ,  2 ,  1 ) 
b a r  (  norm  vecT(  :  ,  normC) ) 

title  ('Norm  col  Tsens,  ABC  Modes  [1:5]'): 

subpl  ot  (3,  2,  3) 

bar(  norm_vecTi  nv(  :  ,  normC) ) 

title  ('Norm  row  TsensINV,  ABC  Modes  [1:5]'); 

subpl  ot  ( 3,  2,  2) 

b  a  r  (  norm  vecTABC(  :  ,  normC)  ) 

title  ('Norm  col  Tsens,  ABC  Modes  [1:10]') 

s  ubpl  ot  (  3,  2,  4) 

b  a  r  ( norm  vecTi  nvABC(  :  ,  normC) ) 

title  ('Norm  row  TsensINV,  ABC  Modes  [1:10]') 

s  ubpl  ot  (  3,  2,  5) 

%  plotting  error  prediction 

%using  [1:5]  modes  of  ABC  system  +[1:5]  modes  of  Base  system; 
baseABCN  =  bar(dv_cal_BasePI  us(:  ,  normC),  .  5,  '  r'  );  hoi  d  on 

%  plotting  error  prediction  using  [1:5]  modes  of  Base  system; 
baseN  =  bar  ( d  v  _  c  a  I  _  A  B  C  (  :  ,  1) ,  .  25,  '  b'  ) ; 


abc_conN  =  i  nt2st  r(  condbasePI  us(  normC)  )  ;  %  cond  no,  for  legend 
abc_conTN  =  spri  ntf(  '  Basel  1:  5]  +ABC[  1:  5]  Cond  =  %s '  ,  abc_conN); 

grid  on 

I  egend(  [  baseABCN,  baseN] ,  BA5ET,  abc_conTN)  ,  hold  on 


%  plotting  actual  error 

if  El  Ibis  -=[  ]  &  mass  Ibis  -=[  ] 

stem(  mass_l  bl  s,  dv  mass.'y',' filled');  hold  on; 
stemf  mass_l  bl  s,  d  v  ~  ma  s  s ,  '  k '  ) 


El  pi  ot  =  El  I  bl  s+10;  hoi  d  on 

stemfEI  Ibis,  dv  EI,'c','filled');hold  on; 

stem(EI_lbis,  d  v  _  E I  ,  '  k 1  ) 


elseif  ma  s  s  _  I  b  I  s  -=[  ]  &  E I  _  I  b  I  s  ==[  ] 

stem(mass_l  bl  s,  dv  mass,'y','filled');hold  on; 
stem(  mass_l  bl  s,  d v ~ ma s s , 1  k '  ) 


else 

stemfEI  Ibis,  d  v  _  E I  ,  '  c'  ,  '  f  i  I  I  e  d '  ) ;  hoi  d  on; 
stem(  El  ~l  bl  s,  d v _ E I  ,  '  k '  ) 


end  %  I  f  E I  _ I  b I  s  ~=[  ]  &  ma s s _ I  b I  s  -=[  ] 


title  (spri  ntf('  Error,  Base  [1:5]  +  ABC  [1:5],  pinned  NODE  #  %d'r  (tt+1))) 


subpl  ot  ( 3,  2,  6) 

%  plotting  error  prediction  using  [1:10]  modes  of  ABC  system; 
baseABCNten  =  bar(  dv_cal  _ABCten(  :  ,  normC)  )  ;  hoi  d  on 
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abc_conNten  =  i  nt  2  s  t  r  ( cond  ABCt  en(  nor  mC )  )  ; 

abc_conTNt  en  =  s  pr  i  nt  f  ( 1  ABt[  1:  10]  Cond  =  %s '  ,  abc_conNten)  ; 

grid  on 

I  egend(  [  baseABCNt  en] ,  abc_conTNt  en) ,  hold  on 


%  plotting  actual  error 

if  El  Ibis  -=[  ]  &  mass  Ibis  -=[  ] 

stem(mass_l  bl  s,  dv  mass.'y'.'filled1);  hold  on; 
stem(mass_lbls,  d v ~ ma s s ,  1  k 1  ) 


El  pi  o t  =  El  I  bl  s+10;  hoi  d  on 

stem(  El  Ibis,  dv  El.'c'.'filled'ljhold  on; 

stem(  El  _l  bl  s,  d v _ E I  ,  1  k 1  ) 


elseif  ma  s  s  _  I  b  I  s  -=[  ]  &  E I  _  I  b  I  s  ==[  ] 

stem(  mass_l  bl  s,  dv  mass, '  y'  ,  'fi  1 1  ed1  );  hoi  d  on; 
stem(  mass_l  bl  s,  d v " ma s s , '  k 1  ) 


else 

stem(  El  _l  bl  s,  d  v  _  E I  ,  1  c'  ,  1  f  i  I  I  ed'  ) ;  hoi  d  on; 
stem(EI"lbis,  d  v  _  E I  ,  1  k '  ) 


end  %  EM  bl  s  -=[  ]  &  mass .1  bl  s  -=[  ] 


title  ( s  pr  i  nt  f  ( 1  Er  r  or ,  ABC  Modes  [1:10],  pinned  NODE  #  %d '  ,  ( 1 1  +1 ) ) ) 

normC  =  normC  +  3;  %  advances  to  the  next  ABC  system 
end  %  1 1  =  1:10  for  plottinggraphes 


END  normRUNthru  crs.m  *********************** 
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p. 


PLOTBEAMMODES  CRS 


********************  plotBea mMo des  crs.m  *********************** 
Calculates  natural  frequencies 
Provided  by  Prof  Gor  di  s 
I  nput  s  needed: 


k  a ,  ma , 

mx ,  kx 

P  r  o  g  r  a  ms 

needed: 

f  Modes 

Out  put  s : 

%  I  a  ma ,  phia,  I  a  mx ,  phix  (without  rigid  body  mo  des) 
%  n  u  m  r  b  m 

%  p  h  i  a  _  p  I  o  t ,  phix  plot 


d  i  s  p  ( '  1  ) ; 

d  i  s  p  ( 1  Calculating  modes  for  each  beam.  ..plot  frequency  comparison1) 


%  Get  modes  of  each  beam: 

[  I  a  ma ,  phi  a  ]  =f  Mo  d  e  s  ( k  a ,  ma )  ; 

[  I  a  mx ,  ph  i  x  ]  =f  Mo  des  ( kx ,  mx ) ; 

%used  to  plot  the  mode  shapes  org  BC  before  ABC 
phia  plot  =  phia; 
phi  x~pl  ot  =  phi  x; 

%  Set  any  rigid  body  mode  freqs  to  zero: 
num_rbm  =  I  ength(fi  nd(l  ama  <  1)); 

spri  ntf('  Number  of  Rigid  Body  Modes  Found:  %2  i  1  ,  num_rbm) 

disp(  '  Removing  rigid  body  mode  frequencies  from  vectors...1) 

I  a  ma  =  I  a  ma  ( f  i  n  d  ( I  a  ma  >  1)5; 

I  a  mx  =  I  a  mx  ( f  i  n  d  ( I  a  mx  >  1 ) ) ; 

o/0  ********************  end  P I  o  t  B  e  a  mMo  d  e  s  crs.m  *********************** 
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Q.  PLOTTINGBARSCRS 


o/0  ********************  pi  ot  t  i  ng  BARS  crs.m  *********************** 

%To  be  used  with  Build2Beams.m  and 

% 

%  This  program  plots  9  graphes  per  figure.  The  first  columns  of  3  graphes 
%  are  the  mode  shapes  of  the  ABC  system  used  in  error  prediction,  The  next 
%  c  o  I  u  mn  of  3  graphes  are  the  error  prediction  using  only  5  modes  of  ABC 
%  system.  The  last  column  of  3  graphes  are  the  error  predictions  using  the 
%  first  5  modes  of  base  system  plus  5  modes  of  the  ABC  system.  The  row 
%  represent  modes  1-5,  middle  row:  modes  6-10,  last  row  :  modes  11-15.  Each 
%  of  the  error  prediction  graphes  also  have  the  base  only  prediction  and 
%  the  actual  error  plotted  for  easy  reference. 

% 

%  Written  by  Constance  R  S  Fernandez,  Spring  2004 

%  Inputs 

% . 

%  c o n d  basePlus 
%  F  OM  ABCSper,  FOM  PLUSper 
%  i  c  n  t  _  o  s  e  t 
%  mo  d  e  s  h  a  p  e 
%  rel  freqERROR 
%  ypos 

%  E I  Ibis,  mass  Ibis 
%  d  v  ’  ma  s  s ,  d  v  _  E  f 

%  Out  put  s 

% . 

%  BASE,  BASET 

%  FOMBASE,  FOMABC,  FOMPLUS 
%  i  n  t  e  r  v  e  I  p 

%  E  R ,  b  a  r  p ,  shape,  error,  a  c  p ,  mo  d  e  p ,  a  p , 

%  model  abel  p,  FOMABCI  abel  p.'FOMPLUSI  abel  p 
%  a  be  con,  abc_conT 
%  ABC"  base 
%  p I  u s _ c o n ,  pi  us_conT 
%  b  a  s  e ,  plus,  E I  _  p  I  o  t 


BASE  =  int2str(cond  bas  eP  I  us  ( 1 ) ) ; 

FOMBASE  =  i  n 1 2 s  t  r  (  F DM  ABC5per (  1)  ) ; 

BASET  =  spri  ntf  (  1  Base'Cond  =  %s ,  FOM  =  %s 1  ,  BASE,  FOMBASE); 


intervelp  =  3; 
mo  d  e  s  h  a  p  e  =  1 ; 

ER  =  1; 

for  barp  =  1:  i  ent  oset 

f  i  gure(  barp+10)  %  figures  11-20 
f  o  r  ma  t  bank 

%s  h  a  p  e  =  [mo  deshape:  mo  deshape +2  0];  %T  h  i  s  number  is  equal  to  number  of  elements, 
shape  =  [  modes  ha  pe:  modes  hape+2  0] ;  %T  h  i  s  number  is  equal  to  number  of  elements, 
error  =  r  o  u  n  d  ( r  e  I  f  r  eq  E  RROR(  E  R;  E  R  +  l  5  )  *  1  0  0  )  /  1  0  0 ; 
a  cp  =  1; 

for  modep  =  1:3  %3  sets  of  modes  per  boundry  condition 

ap  =  [a  cp:  a  cp+4] ;  %modes 
%  REL_errorl  =  int2str(error(a_cp)); 


% 

Errorlabelp  =  spri  ntf(  1  Rel 

error  = 

%s 1  , 

RE  L_ 

er  r  or  1) ; 

% 

REL  e  r  r  o  r  2  =  int2str(error( 

a  cp  +  1) ) 

% 

Errorlabelp  =  spri  ntf('  Rel 

error  = 

'  %s ' 

,  REL 

_  e  r  r  o  r  2 ) 

% 

REL  e  r  r  o  r  3  =  int2str(error( 

a  cp+2) ) 

% 

Errorlabelp  =  spri  ntf('  Rel 

error  = 

'  %s 1 

,  REL 

_  e  r  r  o  r  3 ) 

% 

REL  e  r  r  o  r  4  =  int2str(error( 

a.cp+3) ) 

% 

Errorlabelp  =  spri  ntf('  Rel 

error  = 

’  %s 1 

,  REL 

_  er  r  o r  4 ) 

% 

REL  error5  =  int2str(error( 

a  cp+4)) 

% 

Errorlabelp  =  spri  ntf('  Rel 

error  = 

'  %s ' 

,  REL 

_  e  r  r  o  r  5 ) 

model  abel  p  =  i  n 1 2 s t  r  ( a  c p :  a  cp+4); 

FOMABC  =  i  nt2str(  FOM  ABC5per(  i  ntervel  p  +  modep)  )  ; 
FOMABCI  abel  p  =  spri  ntf('  System  FOM  =  %s 1  ,  FOMABC); 

FOMPLUS  =  int2str(FOM  PLUSper(  i  ntervel  p+modep) ) ; 
FOMPLUSI  abel  p  =  spri  ntf(  1  System  FOM  =  %s 1  ,  FOMPLUS); 
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%  ===============================================================% 

%  =====mode  shape  or  beam  X  and  beam  A  wi  t  h  ABC==================% 

%  ===============================================================% 


fi  gure(  barp+50) 
subpl  o t  ( 3 ,  1,  modep) 
pi  ot  ( y  pos ,  dispA_tot(shape,a 
dispA_tot(shape,a  c  p  + 1 ) , T 
di  spA~t  ot ( shape, a_cp+2) , 
di  spA_tot  ( shape,  a" c p+3) , 
di  spA_tot  ( shape,  a" c p+4) , 
di  spX~tot  ( shape,  a’cp)  , 
di  spX_t  ot  ( shape,  a_cp+l) , 
di  spX_tot  ( shape,  a" c p+2 ) , 
di  spX~t  ot  ( shape,  a_cp+3) , 
di  spX_t  ot  ( shape,  a_cp+4) , 


cp),  ' 

k-  o'  ,  ypos , 

g  -  s 1 , 

ypos,  ... 

b  -  d 1  , 

ypos,  ... 

r  -  x'  , 

ypos,  ... 

m-  * '  , 

ypos,  ... 

r  -  -  o ' 

,  ypos,  .  .  . 

b-  -  s' 

,  ypos,  . . . 

m-  -  d ' 

,  ypos,  . . . 

c  -  -  x ' 

,  ypos  . . . 

k- 

) ,  grid  on. 

I  e  g  e  n  d  (  s  pr  i  nt  f  (  1  Re  I  Freq  Error  =  %d' ,  err  or  (a  c  p ) ) ,  .  .  . 
sprintfl'Rel  Freq  Error  =  %d',  e  r  r  o  r  ( a  c  p  + 1 ) ) ,  .  .  . 

sprintfj'Rel  Freq  Error  =  %d 1  ,  e  r  r  o  r  ( a’  c  p  +2 ) ) ,  .  .  . 

sprintfl'Rel  Freq  Error  =  %d',  error!  a_cp+3)  ) ,  .  .  . 

sprintfl'Rel  Freq  Error  =  %d ' ,  er  r  or  ( a"c  p+4) ) ) ; 

%  I  egendl  s  pr  i  nt  f  ( '  Bm  X,  Md  %d '  ,  a _ c p ' ) ,  s pr  i  nt  f  ( '  Bm  X,  Md  %d '  ,  a  c p  +1 ) ,  . 

%  sprintfl'Bm  X,  Md  %d '  ,  a  cp+2),  s  pr  i  nt  f  ( 1  Bm  X,  Md  %d',a  c  p +3 ) ,  .  .  . 

%  spri  ntf('  Bm  X,  Md  %d '  ,  a ~  c  p  +4 ) ,  spri  ntf(  1  Base,  Md  %d '  ,  a  c  p ) ,  ... 

%  spri  ntf('  Base,  Md  %d '  ,  a"cp+l),  spri  ntf('  Base,  Md  %d '  ,  a"cp+2),  .. 

%  spri  ntf('  Base,  Md  %d '  ,  a"cp+3),  spri  ntf('  Base,  Md  %d '  ,  a"cp+4)) 


t  i  1 1  e(  s pr i  nt f  ( '  Modes  [  %s] '  ,  model  a  be  I  p) ) 
axis  tight 

%  ===============================================================% 

%=====  bar  graphes  of  error  solution  using  only  ABC===  =  ======  =  =  =  % 

%  ===============================================================% 

fi  gure(  barp  +  10)  %  figures  11-20 
subpl  ot(3,  2,  2*modep- 1) 


abc  con  =  int2str(cond  ABC(  i  ntervel  p  +  modep)  )  ; 

abc“conT  =  sprintfl'ABC  Cond  =  %s,  FOM  =  %s 1  ,  a b c _ c o n ,  FOMABC); 

ABC  =  b  a  r  ( d  v  cal  ABC(  :  ,  i  ntervel  p+modep) ,  .  5,  '  r'  ) ;  hold  on 
%A  B  C  =  bar(  d  v  c  a  F  A  B  C  ( : ,  1 : 5 ) ,  .  5 ,  '  r '  ) ;  hold  on  %t  r  y  i  n  g  to  isolate  the  first 
f  i  v  e  mo  d  e  s 

base  =  b  a  r  ( d  v  cal  ABC( : ,  1),  .  25,  '  b'  );  hoi  d  o  f  f  %o  n  %  b  a  s  e  first  5  modes 
%FOMabc  =  barfl,  OF;  hold  off 
grid  on 

%l  egendl  [  ABC,  base,  FOMabc] ,  pi  us  conT,  BASET,  FOMABCI  abel  p) ,  hold  on 
%  g r i  d  on 

I  egendl  [ABC,  base] ,  abc_conT,  BASET) ,  hold  on 
title(sprintf('ABC  only,  [  %s  ] '  ,  model  abel  p)  )  ; 


if  E I  _  I  b  I  s  ~=0  &  mass  Ibis  ~=0 
%  plots  actual  error 

stem!  mass  Ibis,  d  v  mass,  '  y '  ,  'filled');  hold  on; 

stem!  mass“l  bl  s,  d  v  _  ma  s  s , '  k '  ) ,  hold  on; 

El  pi  o  t  =  El  _  I  b  I  s  +1 0 ;  hold  on  %  last  half  of  plot 
s  t  e  m(  E I  pi  o  t ,  dv  EI,'c','filled');hold  on; 
stem(Elplot,  d  v  ’  E I  ,  1  k '  ) ;  hold  on 

%  plots  the  green  triangle  which  indicates  pinned  node 
%pl  ot  ( bar  p+9,  0,  '  gA'  ,  bar  p  +9 ,  0,  1  gh1  ,  bar  p+9,  0,  ’  g*1  ,  .  .  . 

%  barp+11,  0,  1  g'"  ,  barp  +  11,  0,  1  gh'  ,  barp+11,  0,  '  g*'  ) 

plot(((BC(:,3)-l)/2),0,'gA',((BC(:,3)-l)/2),0,'gh',((BC(:,3)- 
1) / 2) ,  0,  '  g* 1  ,  .  .  . 

|(BC(  :  ,  3)-l)/2),  0,  '  gAI  ,  |(BC(:  ,  3)-l)/2),  0,  '  gh'  ,  |(BC(  :  ,  3)-l)/2),  0,  '  g*1 


el  seif  mass  Ibis  ~=0  %&  E I  _  I  b  I  s  =0 
%  p I  o  t  s  actual  error 

stem!  mass  Ibis,  dv  mass,  '  y '  ,  'filled');  hold  on 
stemjmass'lbls,  d  v  ~  ma  s  s ,  '  k '  ) ;  hold  on 
%  plots  the  green  triangle  which  indicates  pinned  node 
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%pl  ot  ( bar  p+9, 0,  1  gA'  ,  bar  p+9,  0,  1  g h '  ,  bar  p+9, 0,  '  g*'  ) 
pi  ot(((BC(  :  ,3)-l)/2),  0,  '  g^1  ,((BC(:  ,  3)-l)/2),0,  1  gh1  ,  ((BC(:  ,  3)- 
1)/ 2) ,  0,  1  g*1  ) 

else 

%  p I  o  t  s  actual  error 

stem)  El  Ibis,  dv  El.'c'.'filled'Jjhold  on; 
stem)  El  Ibis,  d  v  ’  E I  ,  1  k 1  ) ;  hold  on 

%  plots'the  green  triangle  which  indicates  pinned  node 
%pl  ot  ( bar  p+9,  0,  1  gA'  ,  bar  p+9,  0,  1  g h '  ,  bar  p+9,  0,  '  g*'  ) 
pi  ot(((BC(  ;  ,3)-l)/2),  0,  '  g^'  ,((BC(;  ,  3)-l)/2),0,  1  gh1  ,  ((BC(:  ,  3)- 
1)/  2) ,  0,  1  g*1  ) 

end  %  E I  Ibis  ... 


%  ===============================================================% 

%  =========ba r  graphes  of  error  solution  using  ABC  +  ba s e  =  =======% 

%  ===============================================================% 

s  u  b  p  I  o  t  ( 3 ,  2 ,  2  *  mo  d  e  p ) 


plus  con  =  int2str(cond  basePI  us(i  ntervel  p+modep));  %  for  legend 
plus'conT  =  spri  ntf('  Base+ABC  Cond  =  %s  F  O  M  =  %s',  pi  u  s  _  c  o  n ,  F  O  M  P  L  U  S ) ;  %  for 

legend 


plus  =  b  a  r  ( d  v  cal  BasePI  us(:  ,  i  ntervel  p+modep),  .  5,  1  r1  );  hold  on 
base  =  bar(dv'cal'ABC(:,l),.25,'b');  bold  on  %  base  first  5  modes 
%F OMp I  us  =  bar ( 1, D) ;  hold  off 
g  r  i  d  o  n 

legend)  [  pi  us,  base] ,  pi  us  conT,  BASET) ,  hold  on 
%  I  egendj  [  pi  us,  base,  F  OMp  F  us  ]  ,  pi  u  s  _  c  o  n  T ,  BASET,  FOMPLUSI  abel  p) ,  hold  on 

if  E I  _  I  b  I  s  ~=0  &  mass  Ibis  ~=0 
%'  p  I  o  t  s  actual  error 

s  t  e  m(  ma  s  s  Ibis,  dv  ma  s  s ,  1  y '  ,  'filled');  hold  on; 
stemjmass'lbls,  d  v ' ma  s  s ,  1  k '  ) ;  hold  on 
El  pi  o  t  =  E I  _l  b I  s  +1 U ;  hold  on  %  last  half  of  plot 
stem)  El  pi  ot,  dv  El/c'.'filled'jjhold  on; 
stemfElplot,  dv"EI,'k');hold  on 

%  plots  the  green  triangle  which  indicates  pinned  node 
%pl  ot)  bar  p  + 1 ,  0,  1  g A '  ,  bar  p  + 1 ,  0,  1  gh1  ,  bar  p  + 1 ,  0,  1  g*1  ,  .  .  . 

%  barp+11, 0,  1  g^'  ,  barp+11, 0,  1  gh'  ,  barp+11,  0, 1  g*'  ) 
pi  ot  ((( BC(  :  ,  3)  ■  1)  /  2) ,  0,  1  g"  ,  ((  BC(  :  ,  3)- 1)/  2) ,  0,  1  gh1  ,  ( (  BC( ;  ,  3)- 
1)  /  2) ,  0,  1  g*1  ,  .  .  . 

(( BC(  :  ,  3)- 1)  /  2) ,  0,  1  gAI  ,  ( (  BC( :  ,  3)- 1)/  2) ,  0,  1  gh1  ,  ( (  BC(  :  ,  3)-  1)  /  2) ,  0,  '  g*1 


el  sei  f  mass  Ibis  ~=0  %&EI  Ibis  =0 
%  p I  o  t  s  actual  error 

stemfmass  Ibis,  dv  mass,  1  y'  ,  1  fi  I  I  ed1  );  hoi  d  on; 

stem)  mass_l  bl  s,  dv~mass,  '  k'  );  hoi  d  on 

%  plots  the  green  triangle  which  indicates  pinned  node 
%pi  ot)  bar  p  +9 ,  0,  1  gA'  ,  bar  p  +9 ,  0,  1  gh1  ,  bar  p+9,  0,  1  g*1  ) 
plot(((BC(:,3)-l)/2),0,'g/'',((BC(:,3)-l)/2),0,'ghl,((BC(:,3)- 

11/2),  0,  1  g*1  ) 
else 

%  p I  o  t  s  actual  error 

stem)  El  Ibis,  dv  El.'c'.'filled'Jjhold  on; 
stem)  El  ”1  b  I  s ,  d  v  ~  E I  ,  1  k 1  ) ,  hold  on 

%  plots'the  green  triangle  which  indicates  pinned  node 
%  pi  ot)  bar  p  +9 ,  0,  1  gA'  ,  bar  p  +9 ,  0,  1  gh1  ,  bar  p  +9 ,  0,  1  g*1  )%c  hanged  bar  p  + 1 

pi  ot)  (  (  BC(  :  ,  3)-  1)  /  2 ) ,  0,  1  gAI  ,  (  [BC(  :  ,  3)-  1)/  2)  ,  0,  1  gh'  ,  (  (  BC(  :  ,  3)  ■ 

1) / 2) ,  0,  '  g* '  ) 

end  %  i  f  E I _ I  b I  s  ... 

title(sprintf('Base  [1:5]  +  ABC  [  %s  ] 1  ,  model  abel  p)  )  ; 

%  t  i  1 1  e(  s  pr  i  nt  f  ( '  Bas  e  +  ABC,  pinned  at  NODE#  %d 1  ,  b  a  r  p  +1)) 

a  cp=  a _ c p  +5; 
end  %  modep  loop 

mo  deshape  =  mo  deshape  +  11;  %  advances 

interveip  =  intervelp  +3;  %  advances  to  the  next  ABC  system 
ER  =  ER  +  19; 
end  %  b  a  r  p  loop 
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%  f  i  gur  e( bar  p+11) 

%p  I  o  t ( 1 :  2  0 ,  A  B  C 1 ,  bar) 


%  %A B C  Mode  1  Error 
%  subpl  ot(  5,  1,  1) 

%  ABC1  =dv  cal  ABC1 ; 

%  stem)  El  F  b  I  s ,  dv  E I  ,  1  k 1  ) ;  hold  on 

%  b  a  r  (  1:  2D,  ABC1, .  25,  1  r'  )  ; 

%  pi  o t  ( ( (  B C (  :  ,  3 )  - 1 )  /  2 )  ,  0,  1  gA'  ,  ( (  B C (  :  ,  3 )  - 1 )  /  2 ) ,  0,  1  gh'  ,  ( (  B C (  :  ,  3)  - 

1) / 2)  ,  0,  1  g*'  ); 

%  t  i  1 1  e( 1  ABC  Mode  1'  ) 

%  %ABC  Mode  2  Error 
%  subpl  o t  ( 5 ,  1,2) 

% 

%  ABC2  =  dv  cal  ABC2 ; 

%  stem)  El  Ibis,  d  v  E I  ,  1  k 1  ) ;  hold  on 

%  b  a  r  (  1:  2D,  ABC2, .  25, 1  r'  )  ; 

%  pi  ot  (  (( BC(  :  ,  3)  ■  1)  /  2) ,  0,  1  gAI  ,  ( ( BC( :  ,  3)- 1)/  2) ,  0,  '  gh1  ,  ( (  BC(  :  ,  3)- 

1) / 2) ,  0,  '  g* '  )  ; 

%  titlef'ABC  Mode  2') 

% 

%  %A B C  Mode  3  Error 
%  subpl  o t  ( 5 ,  1,3) 

% 

%  ABC3  =  dv  cal  ABC3 ; 

%  stem(  El  Ibis,  d  v  E I  ,  '  k 1  ) ;  hold  on 

%  barf  1:  2D, ABC3, .  25, 1  r' )  ; 

%  pi  ot(((BC(  :  ,  3 )  -  1 )  /  2 )  ,  0,  '  gAI  ,  f(BC(:  ,  3)  -  1 )  /  2 )  ,  0,  '  gh'  ,  f(BC(:  ,  3)  - 

l)/2),0,'g*')  ; 

%  titlef'ABC  Mode  3') 

% 

%  %ABC  Mode  4  Error 
%  s  u  b  p  I  o  t  ( 5 ,  1 , 4 ) 

%  ABC4  =dv  cal  ABC4 ; 

%  stem)  El  Fbls,  dv  E I  ,  '  k 1  ) ;  hold  on 

%  barf  1:  2D, ABC4, .  25, 1  r' )  ; 

%  pi  otfffBCf  :  ,  3 )  -  1 )  /  2 )  ,  0,  '  gAI  ,  ffBCf:  ,  3 )  -  1 )  /  2 )  ,  0,  '  gh'  ,  ffBCf  :  ,  3)  - 

1)/  2)  ,  0,  1  g*'  )  ; 

%  titlef'ABC  Mode  4') 

% 

%  %ABC  Mode  5  Error 
%  s  u  b  p  I  o  t  ( 5 ,  1 ,  5 ) 

%  ABC5  =dv  cal  ABC5 ; 

%  stem)  El  Fbls,  dv  E I  ,  '  k 1  ) ;  hold  on 

%  barf  1:  2D, ABC5, .  25, 1  r' )  ; 

%  pi  otfffBCf  :  ,  3 )  -  1 )  /  2 )  ,  0,  '  gAI  ,  ffBCf:  ,  3)  -  1 )  /  2 )  ,  0,  '  gh'  ,  ffBCf:  ,  3)  - 

l)/2),0,'g*'); 

%  titlef'ABC  Mode  5') 


END  pi  ot  t i  n g BARS  crs.m  *********************** 
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